library(rethinking)
Loading required package: rstan
Loading required package: StanHeaders
Loading required package: ggplot2
Keep up to date with changes at https://www.tidyverse.org/blog/
rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Loading required package: parallel
rethinking (Version 2.11)

Attaching package: ‘rethinking’

The following object is masked from ‘package:stats’:

    rstudent
library(brms)
Loading required package: Rcpp
Loading 'brms' package (version 2.13.5). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: ‘brms’

The following objects are masked from ‘package:rethinking’:

    LOO, stancode, WAIC

The following object is masked from ‘package:rstan’:

    loo

The following object is masked from ‘package:stats’:

    ar
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ tibble  3.0.3     ✓ dplyr   1.0.0
✓ tidyr   1.1.0     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0
✓ purrr   0.3.4     
── Conflicts ───────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x tidyr::extract() masks rstan::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x purrr::map()     masks rethinking::map()
germ <- read_csv("../Dimensions/hydrothermal-all-species/data/light_round1_tall.csv") %>%
  filter(wps == 0) %>%
  select(pops, temps, total_seeds, germ, day, cumulative_germ)
Parsed with column specification:
cols(
  pops = col_character(),
  temps = col_double(),
  wps = col_double(),
  date = col_character(),
  total_seeds = col_double(),
  germ = col_double(),
  start_date = col_date(format = ""),
  census_date = col_date(format = ""),
  day = col_double(),
  cumulative_germ = col_double(),
  cumulative_prop_germ = col_double()
)
germ

what if need one event per row:

one_per_row <- function(df) {
  total_seed <- max(df$total_seeds, sum(df$germ))
  newdata <- tibble(id=1:total_seed, germ=0, day=max(df$day))
  df <- df %>% filter(germ>0)
  count <- 1
  if (nrow(df) > 0) {
    for (i in 1:nrow(df)) { # we look at each row of the df where germination occured
      for (j in 1:df$germ[i]) { # now update the newdata to reflect the germiantion of each seed
        newdata$germ[count] <- 1
        newdata$day[count]=df$day[i]
        count <- count+1 # count keeps track of which individual we are at in the new data
      } # for j
    } # for i
  } # if 
  return(newdata)
}

germone <- germ %>% group_by(pops, temps) %>%
  select(-cumulative_germ) %>% # not needed in this encoding (I think...in any case would need to be recalculated)
  nest() %>%
  mutate(newdata=map(data, one_per_row)) %>%
  select(-data) %>%
  unnest(newdata)

germone

1

Is effect of temp ~linear on cum germ?

germ %>% filter(pops=="STDI", day==28) %>%
  ggplot(aes(x=temps,y=cumulative_germ)) +
  geom_col()

no, so we can’t treat temp as a continuous linear predictor

germ.stdi <- germone %>% filter(pops=="STDI") %>% dplyr::select(-pops)
Adding missing grouping variables: `pops`
germ.stdi

2

censored exponential

rethinking

exponential rate curve, censoring seed that don’t germinate.

d <- list(germ=germ.stdi$germ, 
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

m1.1 <- ulam(
  alist(
    day | germ==1 ~ exponential( lambda),
    day | germ==0 ~ custom(exponential_lccdf( !Y | lambda)),
    lambda <- 1.0 / mu,
    log(mu) <- a[temps],
    a[temps] ~ normal(0,1)),
  data=d,
  chains=4,
  cores = 4
)
starting worker pid=6368 on localhost:11907 at 09:03:25.275
starting worker pid=6382 on localhost:11907 at 09:03:25.610
starting worker pid=6397 on localhost:11907 at 09:03:25.887
starting worker pid=6412 on localhost:11907 at 09:03:26.197

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 8.2e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.82 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 8.6e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.86 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.223272 seconds (Warm-up)
Chain 1:                0.22504 seconds (Sampling)
Chain 1:                0.448312 seconds (Total)
Chain 1: 
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 7.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.278519 seconds (Warm-up)
Chain 2:                0.208053 seconds (Sampling)
Chain 2:                0.486572 seconds (Total)
Chain 2: 
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.274003 seconds (Warm-up)
Chain 3:                0.210096 seconds (Sampling)
Chain 3:                0.484099 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 8.4e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.84 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.260888 seconds (Warm-up)
Chain 4:                0.217804 seconds (Sampling)
Chain 4:                0.478692 seconds (Total)
Chain 4: 
precis(m1.1, depth = 2)

The above represent log(mean time to germination)

exp(5.45)
[1] 232.7582
exp(2.48)
[1] 11.94126

posterior

preddata <- expand_grid(temps=1:8, day=1:28)
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
 $ lambda: num [1:2000, 1:224] 0.00585 0.00443 0.00493 0.0032 0.00844 ...
 $ mu    : num [1:2000, 1:224] 171 226 203 312 119 ...

mu is average days to germ, lambda is rate parameter. neither change over time, of course, so having days doesn’t really make sense.

preddata$mu <- apply(pred$mu,2,mean)
preddata$low <- apply(pred$mu,2,HPDI)[1,]
preddata$high <- apply(pred$mu, 2, HPDI)[2,]

Single temp. Don’t need day posterior

preddata <- expand_grid(temps=3)
pred <- link(m1.1, data = preddata)
str(pred)

how to convert to probs? use pexp.

predprobs <- pexp(1:28,rate=pred$lambda[1])
#plot(x=1:28,y=predprobs, type="l") # crashes on plot

even though it isn’t using day, including day in the prediction data frame will help me keep data in the correct format. Maybe.

preddata <- expand_grid(temps=1:8, day=1:28)
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
 $ lambda: num [1:2000, 1:224] 0.00585 0.00443 0.00493 0.0032 0.00844 ...
 $ mu    : num [1:2000, 1:224] 171 226 203 312 119 ...
predresults <- preddata %>%
  mutate(lambda=as.list(as.data.frame(pred$lambda)))
predresults
predresults <- predresults %>%
  mutate(probs=map2(day, lambda, ~ pexp(.x, .y)),
         mu=map_dbl(probs, mean),
         lower=map_dbl(probs, ~ HPDI(.)[1] %>% unlist()),
         upper=map_dbl(probs, ~ HPDI(.)[2]) %>% unlist())
predresults
predresults %>% select(-lambda, -probs) %>%
  mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
  ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
  geom_line() 

Add realdata:

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

predresults %>% select(-lambda, -probs) %>%
  mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
  ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
  geom_line() +
  geom_point(aes(y=prop_germ), data=stdi.plot)

Poor!

brms censoring model

need to set up indicator for censoring.

germ.stdi <- germ.stdi %>%
  mutate(cens=ifelse(germ==0, "right", "none"),
         tempsc=as.character(temps) %>% str_pad(width=2, pad="0"))
germ.stdi
get_prior(day | cens(cens) ~ 0 + tempsc, family = exponential, data=germ.stdi)
m1.2 <- brm(day | cens(cens) ~ 0 + tempsc,
            family = exponential(),
            set_prior("normal(0,1)", class="b"),
            data = germ.stdi, )
Compiling Stan program...
Start sampling

SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 9.8e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.423825 seconds (Warm-up)
Chain 1:                0.429479 seconds (Sampling)
Chain 1:                0.853304 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 6.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.69 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.426314 seconds (Warm-up)
Chain 2:                0.455533 seconds (Sampling)
Chain 2:                0.881847 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 6.2e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.43039 seconds (Warm-up)
Chain 3:                0.41959 seconds (Sampling)
Chain 3:                0.84998 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL '1de05645edb74de0b2dc64c80aecacfa' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000196 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.96 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.430918 seconds (Warm-up)
Chain 4:                0.426936 seconds (Sampling)
Chain 4:                0.857854 seconds (Total)
Chain 4: 
summary(m1.2)
 Family: exponential 
  Links: mu = log 
Formula: day | cens(cens) ~ 0 + tempsc 
   Data: germ.stdi (Number of observations: 396) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tempsc05     5.44      0.37     4.76     6.22 1.00     6689     2848
tempsc10     4.70      0.28     4.18     5.30 1.00     7327     2780
tempsc15     4.02      0.23     3.60     4.49 1.00     7811     2805
tempsc20     2.48      0.16     2.17     2.81 1.00     7247     2911
tempsc25     3.07      0.17     2.74     3.43 1.00     7072     3064
tempsc30     3.73      0.20     3.36     4.14 1.00     7913     3195
tempsc35     5.21      0.35     4.57     5.94 1.00     6135     3096
tempsc40     5.62      0.40     4.89     6.46 1.00     7016     2782

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
predict(m1.2) %>% head
     Estimate Est.Error     Q2.5     Q97.5
[1,] 249.2430  286.1171 6.667386 1031.5988
[2,] 242.8674  286.0824 5.156902  996.2780
[3,] 252.4897  300.5946 5.439334 1066.8564
[4,] 247.2513  284.7164 5.992272 1018.9538
[5,] 251.3389  286.2790 5.726633 1042.6966
[6,] 243.3522  263.6542 6.032181  953.4948

log(mean time to germination)

predict(m1.2, newdata = expand_grid(tempsc=unique(germ.stdi$tempsc), cens=c("none", "right")))
       Estimate Est.Error      Q2.5      Q97.5
 [1,] 250.54096 285.42410 5.4041250 1048.05543
 [2,] 249.55624 296.61765 5.9321050 1064.49475
 [3,] 115.14309 125.64755 2.5715721  458.46359
 [4,] 114.02507 119.67649 2.7638015  441.18807
 [5,]  57.28506  59.16199 1.4965791  216.83994
 [6,]  57.09426  58.62516 1.4054719  218.42622
 [7,]  12.12453  12.28691 0.3618226   47.03079
 [8,]  12.19735  12.56483 0.3080198   45.58575
 [9,]  21.67803  22.76933 0.5948858   82.15578
[10,]  21.86662  21.88403 0.5093664   78.33396
[11,]  42.78315  43.77602 1.0944283  161.49532
[12,]  43.82032  45.09836 1.1251890  162.88073
[13,] 192.67635 215.60433 3.8413475  773.76129
[14,] 189.78430 217.32587 4.2042633  784.30106
[15,] 300.67681 347.99351 6.5632012 1292.12783
[16,] 293.24047 343.31921 7.0867706 1252.76443
precis(m1.1, depth = 2)

again, these are log(mean time to germination)

lambda should be 1/mu, so

#plot(1:28,pexp(1:28, 1/exp(2.48)), type="l", col="red")
#lines(1:28,pexp(1:28, 1/exp(5.60)), type="l", col="blue")
#crashes
get_prior( day | cens(cens) ~ 0 + tempsc,
           family = Gamma(),
           data = germ.stdi)
m1.3 <- brm(day | cens(cens) ~ 0 + tempsc,
            family = Gamma(),
            set_prior("normal(0,1)", class="b", lb=0),
            data = germ.stdi)
Compiling Stan program...
Start sampling

SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000441 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.41 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 3.83018 seconds (Warm-up)
Chain 1:                2.69174 seconds (Sampling)
Chain 1:                6.52191 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000423 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.23 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 3.7561 seconds (Warm-up)
Chain 2:                2.59898 seconds (Sampling)
Chain 2:                6.35508 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000369 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.69 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 3.65055 seconds (Warm-up)
Chain 3:                2.61099 seconds (Sampling)
Chain 3:                6.26154 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL '77de82ec887792111b6580f2af527efa' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000518 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.18 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 3.6174 seconds (Warm-up)
Chain 4:                2.67537 seconds (Sampling)
Chain 4:                6.29276 seconds (Total)
Chain 4: 
summary(m1.3)
 Family: gamma 
  Links: mu = inverse; shape = identity 
Formula: day | cens(cens) ~ 0 + tempsc 
   Data: germ.stdi (Number of observations: 396) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
tempsc05     0.00      0.00     0.00     0.00 1.00     5074     2617
tempsc10     0.00      0.00     0.00     0.01 1.00     4408     2321
tempsc15     0.01      0.00     0.00     0.02 1.00     4649     2755
tempsc20     0.07      0.02     0.05     0.11 1.00     5340     2742
tempsc25     0.04      0.01     0.02     0.06 1.00     5426     2802
tempsc30     0.02      0.00     0.01     0.03 1.00     5425     2945
tempsc35     0.00      0.00     0.00     0.00 1.00     4444     2688
tempsc40     0.00      0.00     0.00     0.00 1.00     4553     2247

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.68      0.07     0.56     0.82 1.00     3100     2848

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
newdata <- expand_grid(tempsc=unique(germ.stdi$tempsc), cens=c("none", "right"))
cbind(newdata, predict(m1.3, newdata = newdata))

The problem with these models is that they assume that the censored observations are controlled by the same process that is controlling the ones that do germinate. I do not think that is correct, we need a mixture model.

Trying a mixture model for a ZI exponential

look at ulam censored exponential

rethinking::stancode(m1.1)
data{
    int germ[396];
    vector[396] day;
    int temps[396];
}
parameters{
    vector[8] a;
}
model{
    vector[396] lambda;
    vector[396] mu;
    a ~ normal( 0 , 1 );
    for ( i in 1:396 ) {
        mu[i] = a[temps[i]];
        mu[i] = exp(mu[i]);
    }
    for ( i in 1:396 ) {
        lambda[i] = 1/mu[i];
    }
    for ( i in 1:396 ) 
        if ( germ[i] == 0 ) target += exponential_lccdf(day[i] | lambda[i]);
    for ( i in 1:396 ) 
        if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `day`. 

look at brms censored exponential

stancode(m1.2)
// generated with brms 2.13.5
functions {
}
data {
  int<lower=1> N;  // number of observations
  vector[N] Y;  // response variable
  int<lower=-1,upper=2> cens[N];  // indicates censoring
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  vector[K] b;  // population-level effects
}
transformed parameters {
}
model {
  // initialize linear predictor term
  vector[N] mu = X * b;
  for (n in 1:N) {
    // apply the inverse link function
    mu[n] = exp(-(mu[n]));
  }
  // priors including all constants
  target += normal_lpdf(b | 0,1);
  // likelihood including all constants
  if (!prior_only) {
    for (n in 1:N) {
      // special treatment of censored data
      if (cens[n] == 0) {
        target += exponential_lpdf(Y[n] | mu[n]);
      } else if (cens[n] == 1) {
        target += exponential_lccdf(Y[n] | mu[n]);
      } else if (cens[n] == -1) {
        target += exponential_lcdf(Y[n] | mu[n]);
      }
    }
  }
}
generated quantities {
}
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `day`. 
3: Unknown or uninitialised column: `germ`. 
4: Unknown or uninitialised column: `day`. 

Zero inflated, as generated by ulam:

 
There were 18 warnings (use warnings() to see them)
# define parameters
prob_drink <- 0.2 # 20% of days
rate_work <- 1 # average 1 manuscript per day
# sample one year of production 
N <- 365
# simulate days monks drink 
set.seed(365)
drink <- rbinom( N , 1 , prob_drink )
# simulate manuscripts completed
y <- (1-drink)*rpois( N , rate_work )

rethinking::stancode(ulam( alist(
  y|y>0 ~ custom( log1m(p) + poisson_lpmf(y|lambda) ),
  y|y==0 ~ custom( log_mix( p , 0 , poisson_lpmf(0|lambda) ) ),
  logit(p) <- ap,
  log(lambda) <- al,
  ap ~ dnorm(-1.5,1),
  al ~ dnorm(1,0.5)),
  data=list(y=as.integer(y)) , sample=FALSE))
data{
    int y[365];
}
parameters{
    real ap;
    real al;
}
model{
    real p;
    real lambda;
    al ~ normal( 1 , 0.5 );
    ap ~ normal( -1.5 , 1 );
    lambda = al;
    lambda = exp(lambda);
    p = ap;
    p = inv_logit(p);
    for ( i in 1:365 ) 
        if ( y[i] == 0 ) target += log_mix(p, 0, poisson_lpmf(0 | lambda));
    for ( i in 1:365 ) 
        if ( y[i] > 0 ) target += log1m(p) + poisson_lpmf(y[i] | lambda);
}

Zero inflated exponential model, attempt 1

d <- list(N=nrow(germ.stdi),
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `day`. 
3: Unknown or uninitialised column: `germ`. 
4: Unknown or uninitialised column: `day`. 
5: Unknown or uninitialised column: `germ`. 
6: Unknown or uninitialised column: `day`. 
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.4 <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    vector[8] a; // alphas for the linear model
    real ap; // probability of drink
}
model{
    real p;
    vector[N] lambda;
    vector[N] mu;
    a ~ normal( 0 , 1 );
    ap ~ normal( 0, 1 );
    p = ap;
    p = inv_logit(p); // inverse link function for probability drink
    for ( i in 1:N ) {
        mu[i] = a[temps[i]]; // linear model for mu
        mu[i] = exp(mu[i]);  // inverse link function part 1
        lambda[i] = 1/mu[i]; // inverse link function part 2
    }

    for ( i in 1:N ) 
        if ( germ[i] == 0 ) target += log_mix(p, 0, exponential_lccdf(day[i] | lambda[i]));
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) target += log1m(p) +  exponential_lpdf(day[i] | lambda[i]);
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}
"
m1.4 <- stan(model_code=stanmodel1.4, data=d)

SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000112 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.12 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.769353 seconds (Warm-up)
Chain 1:                0.757463 seconds (Sampling)
Chain 1:                1.52682 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000103 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.03 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.855376 seconds (Warm-up)
Chain 2:                0.729975 seconds (Sampling)
Chain 2:                1.58535 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 9.5e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.95 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.798865 seconds (Warm-up)
Chain 3:                0.741906 seconds (Sampling)
Chain 3:                1.54077 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL '883454a76425bb47871cfc2fcde7e236' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 9.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.96 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.801264 seconds (Warm-up)
Chain 4:                0.71224 seconds (Sampling)
Chain 4:                1.5135 seconds (Total)
Chain 4: 
Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-essTail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess
precis(m1.4, depth = 2)
There were 16 warnings (use warnings() to see them)

Model 1.4 fits a single dormancy rate for all, but probably we need to have separate dormancy rate for each temp

take a look at stancode for a binomial model to see what this looks like

db <- list(cumulative_germ=germ$cumulative_germ, temps=as.numeric(as.factor(germ$temps)))
There were 12 warnings (use warnings() to see them)
rethinking::stancode(ulam(alist(
  cumulative_germ ~ dbern(mu),
  mu <- logit(a + b[temps]),
  a ~ dnorm(0,1),
  b[temps] ~ dnorm(0,1)),
  data = db, sample=FALSE, log_lik = TRUE))
data{
    int cumulative_germ[3840];
    int temps[3840];
}
parameters{
    real a;
    vector[8] b;
}
model{
    vector[3840] mu;
    b ~ normal( 0 , 1 );
    a ~ normal( 0 , 1 );
    for ( i in 1:3840 ) {
        mu[i] = logit(a + b[temps[i]]);
    }
    cumulative_germ ~ bernoulli( mu );
}
generated quantities{
    vector[3840] log_lik;
    vector[3840] mu;
    for ( i in 1:3840 ) {
        mu[i] = logit(a + b[temps[i]]);
    }
    for ( i in 1:3840 ) log_lik[i] = bernoulli_lpmf( cumulative_germ[i] | mu[i] );
}

There is no inverse logit because it is built in to the bernoulli function

Try to write a model with different dormancy probs for the different temps.

d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.5 <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    vector[8] a; // alpha for the exponential curve, one for each temp
    vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
    vector[N] p;
    vector[N] lambda;
    vector[N] mu;
    a ~ normal( 0 , 1 );
    ap ~ normal( 0, 1 );
    for (i in 1:N) {
        p[i] = ap[temps[i]];
        p[i] = inv_logit(p[i]); // inverse link function
    }
    for ( i in 1:N ) {
        mu[i] = a[temps[i]];
        mu[i] = exp(mu[i]); // inverse link function
    }
    for ( i in 1:N ) {
        lambda[i] = 1/mu[i];
    }
    for ( i in 1:N ) 
        if ( germ[i] == 0 ) target += log_mix(p[i], 0, exponential_lccdf(day[i] | lambda[i]));
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) target += log1m(p[i]) +  exponential_lpdf(day[i] | lambda[i]);
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}
"
m1.5 <- stan(model_code=stanmodel1.5, data=d)

SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000202 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.02 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.21524 seconds (Warm-up)
Chain 1:                0.957846 seconds (Sampling)
Chain 1:                2.17309 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000121 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.21 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.16162 seconds (Warm-up)
Chain 2:                0.969372 seconds (Sampling)
Chain 2:                2.13099 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000152 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.52 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.18377 seconds (Warm-up)
Chain 3:                0.939891 seconds (Sampling)
Chain 3:                2.12366 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL '9bb5db5579a55e69b320fc828ae7bddf' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.00012 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.17349 seconds (Warm-up)
Chain 4:                0.859005 seconds (Sampling)
Chain 4:                2.03249 seconds (Total)
Chain 4: 
precis(m1.5, depth = 2)
shinystan::launch_shinystan(m1.5)

Launching ShinyStan interface... for large models this  may take some time.

Listening on http://127.0.0.1:3321
post <- as.data.frame(m1.5)
head(post)
post.cor <- cor(post) %>% round(2)
post.cor[abs(post.cor)<.25] <- NA
post.cor
       a[1] a[2] a[3] a[4] a[5] a[6]  a[7] a[8] ap[1] ap[2] ap[3] ap[4] ap[5] ap[6]
a[1]   1.00   NA   NA   NA   NA   NA    NA   NA -0.39    NA    NA    NA    NA    NA
a[2]     NA    1   NA   NA   NA   NA    NA   NA    NA    NA    NA    NA    NA    NA
a[3]     NA   NA    1   NA   NA   NA    NA   NA    NA    NA    NA    NA    NA    NA
a[4]     NA   NA   NA    1   NA   NA    NA   NA    NA    NA    NA    NA    NA    NA
a[5]     NA   NA   NA   NA    1   NA    NA   NA    NA    NA    NA    NA    NA    NA
a[6]     NA   NA   NA   NA   NA    1    NA   NA    NA    NA    NA    NA    NA    NA
a[7]     NA   NA   NA   NA   NA   NA  1.00   NA    NA    NA    NA    NA    NA    NA
a[8]     NA   NA   NA   NA   NA   NA    NA    1    NA    NA    NA    NA    NA    NA
ap[1] -0.39   NA   NA   NA   NA   NA    NA   NA  1.00    NA    NA    NA    NA    NA
ap[2]    NA   NA   NA   NA   NA   NA    NA   NA    NA     1    NA    NA    NA    NA
ap[3]    NA   NA   NA   NA   NA   NA    NA   NA    NA    NA     1    NA    NA    NA
ap[4]    NA   NA   NA   NA   NA   NA    NA   NA    NA    NA    NA     1    NA    NA
ap[5]    NA   NA   NA   NA   NA   NA    NA   NA    NA    NA    NA    NA     1    NA
ap[6]    NA   NA   NA   NA   NA   NA    NA   NA    NA    NA    NA    NA    NA     1
ap[7]    NA   NA   NA   NA   NA   NA -0.32   NA    NA    NA    NA    NA    NA    NA
ap[8]    NA   NA   NA   NA   NA   NA    NA   NA    NA    NA    NA    NA    NA    NA
lp__     NA   NA   NA   NA   NA   NA    NA   NA    NA    NA    NA    NA    NA    NA
      ap[7] ap[8] lp__
a[1]     NA    NA   NA
a[2]     NA    NA   NA
a[3]     NA    NA   NA
a[4]     NA    NA   NA
a[5]     NA    NA   NA
a[6]     NA    NA   NA
a[7]  -0.32    NA   NA
a[8]     NA    NA   NA
ap[1]    NA    NA   NA
ap[2]    NA    NA   NA
ap[3]    NA    NA   NA
ap[4]    NA    NA   NA
ap[5]    NA    NA   NA
ap[6]    NA    NA   NA
ap[7]  1.00    NA   NA
ap[8]    NA     1   NA
lp__     NA    NA    1

now need to make predictions

post_logavgdays <- post %>% select(starts_with("a["
)) # these are the log(mean(avg time to germinate))
mu_lambda <- post_logavgdays %>% 
  summarise_all(mean) %>% # posterior mean
  mutate_all(~ 1/exp(.)) # convert to lambda (inverse link)
mu_lambda
post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(temps=as.factor(unique(germ.stdi$temps)),
                        lambda=t(mu_lambda),
                        p=t(mu_p))

plot, include original data

post_plot <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=pexp(day, rate=lambda) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)

NA

adjust priors

d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.5a <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    vector[8] a; // alpha for the exponential curve, one for each temp
    vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
    vector[N] p;
    vector[N] lambda;
    vector[N] mu;
    a ~ normal( 0 , 2 );
    ap ~ normal( 0, 2 );
    for (i in 1:N) {
        p[i] = ap[temps[i]];
        p[i] = inv_logit(p[i]); // inverse link function
    }
    for ( i in 1:N ) {
        mu[i] = a[temps[i]];
        mu[i] = exp(mu[i]); // inverse link function
    }
    for ( i in 1:N ) {
        lambda[i] = 1/mu[i];
    }
    for ( i in 1:N ) 
        if ( germ[i] == 0 ) target += log_mix(p[i], 0, exponential_lccdf(day[i] | lambda[i]));
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) target += log1m(p[i]) +  exponential_lpdf(day[i] | lambda[i]);
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) day[i] ~ exponential( lambda[i] );
}
"
m1.5a <- stan(model_code=stanmodel1.5a, data=d)

SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00016 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.6 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.2908 seconds (Warm-up)
Chain 1:                0.945991 seconds (Sampling)
Chain 1:                2.23679 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000136 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.36 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1.49546 seconds (Warm-up)
Chain 2:                1.02455 seconds (Sampling)
Chain 2:                2.52001 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000121 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.21 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1.53961 seconds (Warm-up)
Chain 3:                1.05771 seconds (Sampling)
Chain 3:                2.59732 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL '4d1225f5827acf76400029e28fe110a0' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.00012 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1.43043 seconds (Warm-up)
Chain 4:                1.27791 seconds (Sampling)
Chain 4:                2.70834 seconds (Total)
Chain 4: 
There were 1 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problems
precis(m1.5a, depth = 2)
shinystan::launch_shinystan(m1.5a)
post <- as.data.frame(m1.5a)

post_logavgdays <- post %>% select(starts_with("a["
)) # these are the log(mean(avg time to germinate))
mu_lambda <- post_logavgdays %>% 
  summarise_all(mean) %>% # posterior mean
  mutate_all(~ 1/exp(.)) # convert to lambda (inverse link)
mu_lambda

post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p

posterior_coef <- tibble(temps=as.factor(unique(germ.stdi$temps)),
                        lambda=t(mu_lambda),
                        p=t(mu_p))

post_plot <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=pexp(day, rate=lambda) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)

NA

still not great

gamma

Look at gamma stancode from brms:

make_stancode(day ~ tempsc, 
             family="Gamma",
             prior=set_prior("normal(0,1)"),
             data=germ.stdi,
             sample=FALSE)
// generated with brms 2.13.5
functions {
}
data {
  int<lower=1> N;  // number of observations
  vector[N] Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  int Kc = K - 1;
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  }
}
parameters {
  vector[Kc] b;  // population-level effects
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> shape;  // shape parameter
}
transformed parameters {
}
model {
  // initialize linear predictor term
  vector[N] mu = Intercept + Xc * b;
  for (n in 1:N) {
    // apply the inverse link function
    mu[n] = shape * exp(-(mu[n]));
  }
  // priors including all constants
  target += normal_lpdf(b | 0,1);
  target += student_t_lpdf(Intercept | 3, 3.3, 2.5);
  target += gamma_lpdf(shape | 0.01, 0.01);
  // likelihood including all constants
  if (!prior_only) {
    target += gamma_lpdf(Y | shape, mu);
  }
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
}

Try developing a ZI gamma model:

stan uses that gamma(alpha, beta) parameterization where “alpha” is shape" and “beta” is “rate”. So try having the shape be the same for all temps.

d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.7 <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    real<lower=0> shape; // should set lower bound
    vector[8] a; // alpha for the gamma curve, one for each temp
    vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
    vector[N] p;
    vector[N] mu;
    a ~ normal( 0 , .5 ); //narrow priors to overcome divergent transitions
    ap ~ normal( 0, 1.5 ); 
    shape ~ normal(0, .5); //narrow priors to overcome divergent transitions
    for (i in 1:N) {
        p[i] = ap[temps[i]];
        p[i] = inv_logit(p[i]);
    }
    for ( i in 1:N ) {
        mu[i] = a[temps[i]];
         // apply the inverse link function
        mu[i] = shape * exp(-(mu[i]));
    }
    for ( i in 1:N ) 
       if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7 <- stan(model_code=stanmodel1.7, data=d, chains=4, cores=4, control=list(adapt_delta=.99))
starting worker pid=95474 on localhost:11147 at 12:57:13.894
starting worker pid=95488 on localhost:11147 at 12:57:14.340
starting worker pid=95503 on localhost:11147 at 12:57:14.749
starting worker pid=95518 on localhost:11147 at 12:57:15.343

SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001004 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.04 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000521 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 5.21 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.00068 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.8 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '6430b8ea4701020f2821a980ade78a1c' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000521 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.21 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 11.1101 seconds (Warm-up)
Chain 1:                7.85728 seconds (Sampling)
Chain 1:                18.9674 seconds (Total)
Chain 1: 
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 11.5728 seconds (Warm-up)
Chain 2:                7.93563 seconds (Sampling)
Chain 2:                19.5084 seconds (Total)
Chain 2: 
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 12.2292 seconds (Warm-up)
Chain 4:                8.7347 seconds (Sampling)
Chain 4:                20.9639 seconds (Total)
Chain 4: 
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 11.6817 seconds (Warm-up)
Chain 3:                14.9809 seconds (Sampling)
Chain 3:                26.6626 seconds (Total)
Chain 3: 
precis(m1.7, depth = 2)
shinystan::launch_shinystan(m1.7)
post <- as.data.frame(m1.7)
head(post)

now need to make predictions

post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>% 
  mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
  summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)
mu_shape
[1] 2.350936
post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p

plot predictions and original data

posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
                        rate=t(mu_rate),
                        p=t(mu_p))

post_plot1.7 <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot1.7 %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)

NA

playing with priors:

Make shape exponential

d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.7a <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    real<lower=0> shape; // should set lower bound
    vector[8] a; // alpha for the gamma curve, one for each temp
    vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
    vector[N] p;
    vector[N] mu;
    a ~ normal( 0 , .5 ); //narrow priors to overcome divergent transitions
    ap ~ normal( 0, 1.5 ); 
    shape ~ exponential(.5); 
    for (i in 1:N) {
        p[i] = ap[temps[i]];
        p[i] = inv_logit(p[i]);
    }
    for ( i in 1:N ) {
        mu[i] = a[temps[i]];
         // apply the inverse link function
        mu[i] = shape * exp(-(mu[i]));
    }
    for ( i in 1:N ) 
       if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7a <- stan(model_code=stanmodel1.7a, data=d, chains=4, cores=4, control=list(adapt_delta=.99))
starting worker pid=8458 on localhost:11686 at 09:35:49.473
starting worker pid=8472 on localhost:11686 at 09:35:49.776
starting worker pid=8487 on localhost:11686 at 09:35:50.142
starting worker pid=8502 on localhost:11686 at 09:35:50.398

SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Gradient evaluated at the initial value is not finite.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Gradient evaluation took 0.000855 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.55 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000976 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.76 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000618 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '3c231baf7e3835827dc1a0e1dda6a809' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000481 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.81 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 10.5168 seconds (Warm-up)
Chain 1:                7.28403 seconds (Sampling)
Chain 1:                17.8008 seconds (Total)
Chain 1: 
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 11.3064 seconds (Warm-up)
Chain 3:                7.1714 seconds (Sampling)
Chain 3:                18.4778 seconds (Total)
Chain 3: 
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 14.172 seconds (Warm-up)
Chain 4:                7.23646 seconds (Sampling)
Chain 4:                21.4084 seconds (Total)
Chain 4: 
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 10.5926 seconds (Warm-up)
Chain 2:                14.0434 seconds (Sampling)
Chain 2:                24.636 seconds (Total)
Chain 2: 
precis(m1.7a, depth = 2)
post <- as.data.frame(m1.7a)

post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>% 
  mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
  summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)

post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p

posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
                        rate=t(mu_rate),
                        p=t(mu_p))

post_plot1.7a <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot1.7a %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)

NA

Shape parameter is larger than in 1.7. The fit is slightly better at least for 20 and 25.

make a comparison plot

post_plot_compare <- bind_rows(old=post_plot1.7, new=post_plot1.7a, .id = "fit") 
post_plot_compare %>%
  ggplot(aes(x=day,y=prop_germ,color=temps,group=str_c(temps,fit),linetype=fit)) +
  geom_line() +
  geom_point(data=stdi.plot, aes(x=day,y=prop_germ,color=temps), inherit.aes = FALSE)

NA

Just a hair better…

continue to work with priors

d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.7b <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    real<lower=0> shape; // should set lower bound
    vector[8] a; // alpha for the gamma curve, one for each temp
    vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
    vector[N] p;
    vector[N] mu;
    a ~ normal( 0 , .75 ); // divergent with 1
    ap ~ normal( 0, 1.5 ); 
    shape ~ exponential(1.5); // divergent if this is .5 and a is 0, 1
    for (i in 1:N) {
        p[i] = ap[temps[i]];
        p[i] = inv_logit(p[i]);
    }
    for ( i in 1:N ) {
        mu[i] = a[temps[i]];
         // apply the inverse link function
        mu[i] = shape * exp(-(mu[i]));
    }
    for ( i in 1:N ) 
       if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7b <- stan(model_code=stanmodel1.7b, data=d, chains=4, cores=4, control=list(adapt_delta=.99))
starting worker pid=8771 on localhost:11686 at 09:42:59.153
starting worker pid=8785 on localhost:11686 at 09:42:59.459
starting worker pid=8800 on localhost:11686 at 09:42:59.791
starting worker pid=8815 on localhost:11686 at 09:43:00.136

SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000634 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.34 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000827 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 8.27 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000618 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 6.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'a70b71c705816de2eb4120643abe8088' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.00093 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 9.3 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 12.1323 seconds (Warm-up)
Chain 1:                7.61628 seconds (Sampling)
Chain 1:                19.7486 seconds (Total)
Chain 1: 
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 15.0125 seconds (Warm-up)
Chain 4:                7.62409 seconds (Sampling)
Chain 4:                22.6365 seconds (Total)
Chain 4: 
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 12.2946 seconds (Warm-up)
Chain 3:                12.3514 seconds (Sampling)
Chain 3:                24.6461 seconds (Total)
Chain 3: 
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 12.6255 seconds (Warm-up)
Chain 2:                14.3954 seconds (Sampling)
Chain 2:                27.0209 seconds (Total)
Chain 2: 
There were 39 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problems
precis(m1.7b, depth = 2)
post <- as.data.frame(m1.7b)

post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>% 
  mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
  summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)

post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p

posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
                        rate=t(mu_rate),
                        p=t(mu_p))

post_plot <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)

NA

no improvement over 1.7a

brms hurdle gamma

get_prior(bf(day ~ 0 + tempsc, hu ~ 0 + tempsc),
          family = hurdle_gamma,
          data=germ.stdi)
make_stancode(bf(day ~ 0 + tempsc, hu ~ 0 + tempsc),
          family = hurdle_gamma,
          prior=set_prior("normal(0,1)"),
          data=germ.stdi)
m1.8 <- brm(bf(day ~ 0 + tempsc, hu ~ 0 + tempsc),
          family = hurdle_gamma,
          prior=c(set_prior("normal(0,1)"),
                  set_prior("normal(0,1)", dpar="hu")), # must specify this!
          data=germ.stdi, cores = 4)
m1.8
post <- as.data.frame(m1.8)
head(post)

weird fit!

now need to make predictions

post_logavgdays <- post %>% select(shape, starts_with("b_temp"
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>% 
  mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
  summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)
post_logitp <- post %>% select(starts_with("b_hu")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
                        rate=t(mu_rate),
                        p=t(mu_p))

plot predictions and original data



stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)
  

allow both shape and rate of ZI gamma to be predicted by temp

Don’t run. Doesn’t sample well because shape and rate are correlated.

d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.9 <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    vector[8] as; // alpha for gama shape
    vector[8] ar; // alpha for the gamma rate, one for each temp
    vector[8] ap; // alpha for the proportion dormant, one for each temp
}
model{
    vector[N] p;
    vector[N] shape;
    vector[N] rate;
    as ~ normal( 0, 1 );
    ar ~ normal( 0, 1 );
    ap ~ normal( 0, 1.5 ); 
    
    for (i in 1:N) {
        p[i] = ap[temps[i]];
        p[i] = inv_logit(p[i]);
    }
    for ( i in 1:N ) {
        shape[i] = as[temps[i]];
        rate[i] = ar[temps[i]];
         // apply the inverse link functions
        shape[i] = exp(shape[i]);
        rate[i] = exp(rate[i]);
    }
    for ( i in 1:N ) 
       if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape[i], rate[i]));
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape[i], rate[i]);
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) day[i] ~ gamma( shape[i], rate [i] );
}
generated quantities {

}
"
m1.9 <- stan(model_code=stanmodel1.9, data=d, chains=4, cores=4, control=list(adapt_delta=.99))
precis(m1.9, depth = 2)
post <- as.data.frame(m1.9)

post_logavgdays <- post %>% select(shape, starts_with("a["
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>% 
  mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
  summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)

post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p

posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
                        rate=t(mu_rate),
                        p=t(mu_p))

post_plot <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)
  

log normal

make_stancode(day ~ tempsc, family = lognormal(), data=germ.stdi)
// generated with brms 2.13.5
functions {
}
data {
  int<lower=1> N;  // number of observations
  vector[N] Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  int Kc = K - 1;
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  }
}
parameters {
  vector[Kc] b;  // population-level effects
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> sigma;  // residual SD
}
transformed parameters {
}
model {
  // initialize linear predictor term
  vector[N] mu = Intercept + Xc * b;
  // priors including all constants
  target += student_t_lpdf(Intercept | 3, 3.3, 2.5);
  target += student_t_lpdf(sigma | 3, 0, 2.5)
    - 1 * student_t_lccdf(0 | 3, 0, 2.5);
  // likelihood including all constants
  if (!prior_only) {
    target += lognormal_lpdf(Y | mu, sigma);
  }
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
}
d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.10 <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    vector[8] a; // alpha for the lognormal curve, one for each temp
    vector[8] ap; // alpha for the proportion dormant, one for each temp
    real<lower=0> sigma;
}
model{
    vector[N] p;
    vector[N] mu;
    a ~ normal( 0 , 2 );
    sigma ~ exponential (1);
    ap ~ normal( 0, 2 );
    for (i in 1:N) {
        p[i] = ap[temps[i]];
        p[i] = inv_logit(p[i]); // inverse link function
    }
    for ( i in 1:N ) {
        mu[i] = a[temps[i]];
    }
    for ( i in 1:N ) 
        if ( germ[i] == 0 ) target += log_mix(p[i], 0, lognormal_lccdf(day[i] | mu[i], sigma));
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) target += log1m(p[i]) +  lognormal_lpdf(day[i] | mu[i], sigma);
    for ( i in 1:N ) 
        if ( germ[i] == 1 ) day[i] ~ lognormal( mu[i], sigma );
}
"
m1.10 <- stan(model_code=stanmodel1.10, data=d)

SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000125 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.25 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 3.14968 seconds (Warm-up)
Chain 1:                2.42543 seconds (Sampling)
Chain 1:                5.57511 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000141 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.41 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 2.67482 seconds (Warm-up)
Chain 2:                1.71738 seconds (Sampling)
Chain 2:                4.3922 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000107 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.07 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 2.77018 seconds (Warm-up)
Chain 3:                1.81761 seconds (Sampling)
Chain 3:                4.58779 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'aef830d7cd66c8976ff411b35f013384' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000139 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.39 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 3.14339 seconds (Warm-up)
Chain 4:                3.44874 seconds (Sampling)
Chain 4:                6.59213 seconds (Total)
Chain 4: 
There were 35 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problems
precis(m1.10, depth = 2)
shinystan::launch_shinystan(m1.10)
post <- as.data.frame(m1.10)
head(post)

now need to make predictions

post_logavgdays <- post %>% select(starts_with("a["
)) # these are the log(mean(avg time to germinate))
mu <- post_logavgdays %>% 
  summarise_all(mean) 

mu_sigma <- mean(post$sigma)
post_logitp <- post %>% select(starts_with("ap")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p
posterior_coef <- tibble(temps=as.factor(unique(germ.stdi$temps)),
                        mu=t(mu),
                        p=t(mu_p),
                        sigma=mu_sigma)

plot, include original data

post_plot <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=plnorm(day, meanlog=mu, sdlog=sigma) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)

NA

week2

mv norm for rate and p?

Look at Stancode from chapter 14 mvnorm

rethinking::stancode(ulam(
    alist(
        wait ~ normal( mu , sigma ),
        mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
        c(a_cafe,b_cafe)[cafe] ~ multi_normal( c(a,b) , Rho , sigma_cafe ),
        a ~ normal(5,2),
        b ~ normal(-1,0.5),
        sigma_cafe ~ exponential(1),
        sigma ~ exponential(1),
        Rho ~ lkj_corr(2)
    ) , data=d , sample = FALSE ))
data{
    vector[200] wait;
    int afternoon[200];
    int cafe[200];
}
parameters{
    vector[20] b_cafe;
    vector[20] a_cafe;
    real a;
    real b;
    vector<lower=0>[2] sigma_cafe;
    real<lower=0> sigma;
    corr_matrix[2] Rho;
}
model{
    vector[200] mu;
    Rho ~ lkj_corr( 2 );
    sigma ~ exponential( 1 );
    sigma_cafe ~ exponential( 1 );
    b ~ normal( -1 , 0.5 );
    a ~ normal( 5 , 2 );
    {
    vector[2] YY[20];
    vector[2] MU;
    MU = [ a , b ]';
    for ( j in 1:20 ) YY[j] = [ a_cafe[j] , b_cafe[j] ]';
    YY ~ multi_normal( MU , quad_form_diag(Rho , sigma_cafe) );
    }
    for ( i in 1:200 ) {
        mu[i] = a_cafe[cafe[i]] + b_cafe[cafe[i]] * afternoon[i];
    }
    wait ~ normal( mu , sigma );
}

Make shape exponential

d <- list(N=nrow(germ.stdi),
          germ=germ.stdi$germ,
          temps=as.numeric(as.factor(germ.stdi$temps)),
          day=germ.stdi$day)

stanmodel1.7mv <-
  "
data{
    int<lower=1> N;  // number of observations
    int germ[N];
    vector[N] day;
    int temps[N];
}
parameters{
    real<lower=0> shape; // should set lower bound
    vector[8] a_temp; // alpha for the gamma curve, one for each temp
    vector[8] ap_temp; // alpha for the proportion dormant, one for each temp
    real a;
    real ap;
    vector<lower=0>[2] sigma_temp;
    corr_matrix[2] Rho;
}
model{
    vector[N] p;
    vector[N] mu;
    Rho ~ lkj_corr(2);
    sigma_temp ~ exponential(1);
    
    a ~ normal( 0 , .5 ); //narrow priors to overcome divergent transitions
    ap ~ normal( 0, 1.5 ); 
    shape ~ exponential(.5); 
{
    vector[2] YY[8];
    vector[2] MU;
    MU = [a, ap]';
    for (j in 1:8) YY[j] = [a_temp[j], ap_temp[j]]';
    YY ~ multi_normal(MU, quad_form_diag(Rho, sigma_temp));
}
    for (i in 1:N) {
        p[i] = ap_temp[temps[i]];
        p[i] = inv_logit(p[i]);
    }
    for ( i in 1:N ) {
        mu[i] = a_temp[temps[i]];
         // apply the inverse link function
        mu[i] = shape * exp(-(mu[i]));
    }
    for ( i in 1:N ) 
       if ( germ[i] == 0 ) target += log_mix(p[i], 0, gamma_lccdf(day[i] | shape, mu[i]));
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) target += log1m(p[i]) + gamma_lpdf(day[i] | shape, mu[i]);
    for ( i in 1:N ) 
       if ( germ[i] == 1 ) day[i] ~ gamma( shape, mu[i] );
}
"
m1.7mv <- stan(model_code=stanmodel1.7mv, data=d, chains=4, cores=4, control=list(adapt_delta=.99))
starting worker pid=95107 on localhost:11147 at 12:45:01.786
starting worker pid=95122 on localhost:11147 at 12:45:02.133
starting worker pid=95137 on localhost:11147 at 12:45:02.451
starting worker pid=95152 on localhost:11147 at 12:45:02.774

SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000776 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.76 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.001002 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 10.02 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000718 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL '3c88bdab9f8e99a16782291c3b0f1b59' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000593 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 5.93 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 27.0793 seconds (Warm-up)
Chain 1:                26.5605 seconds (Sampling)
Chain 1:                53.6397 seconds (Total)
Chain 1: 
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 27.0663 seconds (Warm-up)
Chain 3:                27.5972 seconds (Sampling)
Chain 3:                54.6635 seconds (Total)
Chain 3: 
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 27.7892 seconds (Warm-up)
Chain 2:                28.2471 seconds (Sampling)
Chain 2:                56.0363 seconds (Total)
Chain 2: 
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 34.5539 seconds (Warm-up)
Chain 4:                22.4221 seconds (Sampling)
Chain 4:                56.976 seconds (Total)
Chain 4: 
There were 100 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Examine the pairs() plot to diagnose sampling problems

100 divergent transitions…

precis(m1.7mv, depth = 3)
post <- as.data.frame(m1.7mv)

post_logavgdays <- post %>% select(shape, starts_with("a_"
)) # these are the shape and the log(mean(avg time to germinate))
mu_rate <- post_logavgdays %>% 
  mutate(across(-shape, ~ shape*exp(-(.))) ) %>%
  summarize(across(everything(), mean)) %>% select(-shape) # posterior mean
mu_rate
mu_shape <- mean(post$shape)

post_logitp <- post %>% select(starts_with("ap_")) #logit p dormant
mu_p <- post_logitp %>%
  summarize_all(mean) %>%
  mutate_all(inv_logit)
mu_p

posterior_coef <- tibble(shape=mu_shape, temps=as.factor(unique(germ.stdi$temps)),
                        rate=t(mu_rate),
                        p=t(mu_p))

post_plot1.7mv <- expand_grid(posterior_coef, day=1:28) %>%
  mutate(prop_germ=pgamma(day, shape=mu_shape, rate=rate) * (1-p))

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop_germ=cumulative_germ/total_seeds)

post_plot1.7mv %>% 
  ggplot(aes(x=day,y=prop_germ,color=temps,group=temps)) +
  geom_line() +
  geom_point(data=stdi.plot)

NA

make a comparison plot

post_plot_compare <- bind_rows(old=post_plot1.7, new=post_plot1.7mv, .id = "fit") 
post_plot_compare %>%
  ggplot(aes(x=day,y=prop_germ,color=temps,group=str_c(temps,fit),linetype=fit)) +
  geom_line() +
  geom_point(data=stdi.plot, aes(x=day,y=prop_germ,color=temps), inherit.aes = FALSE)

NA

overall not as good

extend to multiple species/pops

modify model 1.7a

length(unique(germone$pops))
[1] 24

how can we switch this around to use day as a predictor?

LS0tCnRpdGxlOiAiZ2VybWluYXRpb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyfQpsaWJyYXJ5KHJldGhpbmtpbmcpCmxpYnJhcnkoYnJtcykKbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKYGBge3J9Cmdlcm0gPC0gcmVhZF9jc3YoIi4uL0RpbWVuc2lvbnMvaHlkcm90aGVybWFsLWFsbC1zcGVjaWVzL2RhdGEvbGlnaHRfcm91bmQxX3RhbGwuY3N2IikgJT4lCiAgZmlsdGVyKHdwcyA9PSAwKSAlPiUKICBzZWxlY3QocG9wcywgdGVtcHMsIHRvdGFsX3NlZWRzLCBnZXJtLCBkYXksIGN1bXVsYXRpdmVfZ2VybSkKZ2VybQpgYGAKCndoYXQgaWYgbmVlZCBvbmUgZXZlbnQgcGVyIHJvdzoKCmBgYHtyfQpvbmVfcGVyX3JvdyA8LSBmdW5jdGlvbihkZikgewogIHRvdGFsX3NlZWQgPC0gbWF4KGRmJHRvdGFsX3NlZWRzLCBzdW0oZGYkZ2VybSkpCiAgbmV3ZGF0YSA8LSB0aWJibGUoaWQ9MTp0b3RhbF9zZWVkLCBnZXJtPTAsIGRheT1tYXgoZGYkZGF5KSkKICBkZiA8LSBkZiAlPiUgZmlsdGVyKGdlcm0+MCkKICBjb3VudCA8LSAxCiAgaWYgKG5yb3coZGYpID4gMCkgewogICAgZm9yIChpIGluIDE6bnJvdyhkZikpIHsgIyB3ZSBsb29rIGF0IGVhY2ggcm93IG9mIHRoZSBkZiB3aGVyZSBnZXJtaW5hdGlvbiBvY2N1cmVkCiAgICAgIGZvciAoaiBpbiAxOmRmJGdlcm1baV0pIHsgIyBub3cgdXBkYXRlIHRoZSBuZXdkYXRhIHRvIHJlZmxlY3QgdGhlIGdlcm1pYW50aW9uIG9mIGVhY2ggc2VlZAogICAgICAgIG5ld2RhdGEkZ2VybVtjb3VudF0gPC0gMQogICAgICAgIG5ld2RhdGEkZGF5W2NvdW50XT1kZiRkYXlbaV0KICAgICAgICBjb3VudCA8LSBjb3VudCsxICMgY291bnQga2VlcHMgdHJhY2sgb2Ygd2hpY2ggaW5kaXZpZHVhbCB3ZSBhcmUgYXQgaW4gdGhlIG5ldyBkYXRhCiAgICAgIH0gIyBmb3IgagogICAgfSAjIGZvciBpCiAgfSAjIGlmIAogIHJldHVybihuZXdkYXRhKQp9CgpnZXJtb25lIDwtIGdlcm0gJT4lIGdyb3VwX2J5KHBvcHMsIHRlbXBzKSAlPiUKICBzZWxlY3QoLWN1bXVsYXRpdmVfZ2VybSkgJT4lICMgbm90IG5lZWRlZCBpbiB0aGlzIGVuY29kaW5nIChJIHRoaW5rLi4uaW4gYW55IGNhc2Ugd291bGQgbmVlZCB0byBiZSByZWNhbGN1bGF0ZWQpCiAgbmVzdCgpICU+JQogIG11dGF0ZShuZXdkYXRhPW1hcChkYXRhLCBvbmVfcGVyX3JvdykpICU+JQogIHNlbGVjdCgtZGF0YSkgJT4lCiAgdW5uZXN0KG5ld2RhdGEpCgpnZXJtb25lCmBgYAoKIyMgMQoKSXMgZWZmZWN0IG9mIHRlbXAgfmxpbmVhciBvbiBjdW0gZ2VybT8KCmBgYHtyfQpnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiLCBkYXk9PTI4KSAlPiUKICBnZ3Bsb3QoYWVzKHg9dGVtcHMseT1jdW11bGF0aXZlX2dlcm0pKSArCiAgZ2VvbV9jb2woKQpgYGAKbm8sIHNvIHdlIGNhbid0IHRyZWF0IHRlbXAgYXMgYSBjb250aW51b3VzIGxpbmVhciBwcmVkaWN0b3IKCmBgYHtyfQpnZXJtLnN0ZGkgPC0gZ2VybW9uZSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIGRwbHlyOjpzZWxlY3QoLXBvcHMpCmdlcm0uc3RkaQpgYGAKCiMjIDIKCiMjIyBjZW5zb3JlZCBleHBvbmVudGlhbAoKIyMjIyByZXRoaW5raW5nCgpleHBvbmVudGlhbCByYXRlIGN1cnZlLCBjZW5zb3Jpbmcgc2VlZCB0aGF0IGRvbid0IGdlcm1pbmF0ZS4KYGBge3J9CmQgPC0gbGlzdChnZXJtPWdlcm0uc3RkaSRnZXJtLCAKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgptMS4xIDwtIHVsYW0oCiAgYWxpc3QoCiAgICBkYXkgfCBnZXJtPT0xIH4gZXhwb25lbnRpYWwoIGxhbWJkYSksCiAgICBkYXkgfCBnZXJtPT0wIH4gY3VzdG9tKGV4cG9uZW50aWFsX2xjY2RmKCAhWSB8IGxhbWJkYSkpLAogICAgbGFtYmRhIDwtIDEuMCAvIG11LAogICAgbG9nKG11KSA8LSBhW3RlbXBzXSwKICAgIGFbdGVtcHNdIH4gbm9ybWFsKDAsMSkpLAogIGRhdGE9ZCwKICBjaGFpbnM9NCwKICBjb3JlcyA9IDQKKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuMSwgZGVwdGggPSAyKQpgYGAKVGhlIGFib3ZlIHJlcHJlc2VudCBsb2cobWVhbiB0aW1lIHRvIGdlcm1pbmF0aW9uKQoKYGBge3J9CmV4cCg1LjQ1KQpleHAoMi40OCkKYGBgCgoKcG9zdGVyaW9yCmBgYHtyfQpwcmVkZGF0YSA8LSBleHBhbmRfZ3JpZCh0ZW1wcz0xOjgsIGRheT0xOjI4KQpwcmVkIDwtIGxpbmsobTEuMSwgZGF0YSA9IHByZWRkYXRhKQpzdHIocHJlZCkKYGBgCgptdSBpcyBhdmVyYWdlIGRheXMgdG8gZ2VybSwgbGFtYmRhIGlzIHJhdGUgcGFyYW1ldGVyLiAgbmVpdGhlciBjaGFuZ2Ugb3ZlciB0aW1lLCBvZiBjb3Vyc2UsIHNvIGhhdmluZyBkYXlzIGRvZXNuJ3QgcmVhbGx5IG1ha2Ugc2Vuc2UuCgpgYGB7cn0KcHJlZGRhdGEkbXUgPC0gYXBwbHkocHJlZCRtdSwyLG1lYW4pCnByZWRkYXRhJGxvdyA8LSBhcHBseShwcmVkJG11LDIsSFBESSlbMSxdCnByZWRkYXRhJGhpZ2ggPC0gYXBwbHkocHJlZCRtdSwgMiwgSFBESSlbMixdCmBgYAoKU2luZ2xlIHRlbXAuICBEb24ndCBuZWVkIGRheQpwb3N0ZXJpb3IKYGBge3J9CnByZWRkYXRhIDwtIGV4cGFuZF9ncmlkKHRlbXBzPTMpCnByZWQgPC0gbGluayhtMS4xLCBkYXRhID0gcHJlZGRhdGEpCnN0cihwcmVkKQpgYGAKCmhvdyB0byBjb252ZXJ0IHRvIHByb2JzPyB1c2UgcGV4cC4KCmBgYHtyfQpwcmVkcHJvYnMgPC0gcGV4cCgxOjI4LHJhdGU9cHJlZCRsYW1iZGFbMV0pCmBgYAoKCmBgYHtyfQojcGxvdCh4PTE6MjgseT1wcmVkcHJvYnMsIHR5cGU9ImwiKSAjIGNyYXNoZXMgb24gcGxvdApgYGAKCmV2ZW4gdGhvdWdoIGl0IGlzbid0IHVzaW5nIGRheSwgaW5jbHVkaW5nIGRheSBpbiB0aGUgcHJlZGljdGlvbiBkYXRhIGZyYW1lIHdpbGwgaGVscCBtZSBrZWVwIGRhdGEgaW4gdGhlIGNvcnJlY3QgZm9ybWF0LiAgTWF5YmUuCgpgYGB7cn0KcHJlZGRhdGEgPC0gZXhwYW5kX2dyaWQodGVtcHM9MTo4LCBkYXk9MToyOCkKcHJlZCA8LSBsaW5rKG0xLjEsIGRhdGEgPSBwcmVkZGF0YSkKc3RyKHByZWQpCmBgYAoKYGBge3J9CnByZWRyZXN1bHRzIDwtIHByZWRkYXRhICU+JQogIG11dGF0ZShsYW1iZGE9YXMubGlzdChhcy5kYXRhLmZyYW1lKHByZWQkbGFtYmRhKSkpCnByZWRyZXN1bHRzCmBgYAoKYGBge3J9CnByZWRyZXN1bHRzIDwtIHByZWRyZXN1bHRzICU+JQogIG11dGF0ZShwcm9icz1tYXAyKGRheSwgbGFtYmRhLCB+IHBleHAoLngsIC55KSksCiAgICAgICAgIG11PW1hcF9kYmwocHJvYnMsIG1lYW4pLAogICAgICAgICBsb3dlcj1tYXBfZGJsKHByb2JzLCB+IEhQREkoLilbMV0gJT4lIHVubGlzdCgpKSwKICAgICAgICAgdXBwZXI9bWFwX2RibChwcm9icywgfiBIUERJKC4pWzJdKSAlPiUgdW5saXN0KCkpCnByZWRyZXN1bHRzCmBgYAoKCmBgYHtyfQpwcmVkcmVzdWx0cyAlPiUgc2VsZWN0KC1sYW1iZGEsIC1wcm9icykgJT4lCiAgbXV0YXRlKHRlbXBzPWZhY3Rvcih0ZW1wcywgbGFiZWxzPWFzLmNoYXJhY3Rlcihzb3J0KHVuaXF1ZShnZXJtLnN0ZGkkdGVtcHMpKSkpKSAlPiUKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9bXUsY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgCmBgYAoKQWRkIHJlYWxkYXRhOgoKYGBge3J9CnN0ZGkucGxvdCA8LSBnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiKSAlPiUgCiAgc2VsZWN0KGRheSwgdGVtcHMsIGN1bXVsYXRpdmVfZ2VybSwgdG90YWxfc2VlZHMpICU+JQogIG11dGF0ZSh0ZW1wcz1hcy5mYWN0b3IodGVtcHMpLAogICAgICAgICBwcm9wX2dlcm09Y3VtdWxhdGl2ZV9nZXJtL3RvdGFsX3NlZWRzKQoKcHJlZHJlc3VsdHMgJT4lIHNlbGVjdCgtbGFtYmRhLCAtcHJvYnMpICU+JQogIG11dGF0ZSh0ZW1wcz1mYWN0b3IodGVtcHMsIGxhYmVscz1hcy5jaGFyYWN0ZXIoc29ydCh1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSkpKSkgJT4lCiAgZ2dwbG90KGFlcyh4PWRheSx5PW11LGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGFlcyh5PXByb3BfZ2VybSksIGRhdGE9c3RkaS5wbG90KQpgYGAKClBvb3IhCgojIyMgYnJtcyBjZW5zb3JpbmcgbW9kZWwKCm5lZWQgdG8gc2V0IHVwIGluZGljYXRvciBmb3IgY2Vuc29yaW5nLgpgYGB7cn0KZ2VybS5zdGRpIDwtIGdlcm0uc3RkaSAlPiUKICBtdXRhdGUoY2Vucz1pZmVsc2UoZ2VybT09MCwgInJpZ2h0IiwgIm5vbmUiKSwKICAgICAgICAgdGVtcHNjPWFzLmNoYXJhY3Rlcih0ZW1wcykgJT4lIHN0cl9wYWQod2lkdGg9MiwgcGFkPSIwIikpCmdlcm0uc3RkaQpgYGAKCmBgYHtyfQpnZXRfcHJpb3IoZGF5IHwgY2VucyhjZW5zKSB+IDAgKyB0ZW1wc2MsIGZhbWlseSA9IGV4cG9uZW50aWFsLCBkYXRhPWdlcm0uc3RkaSkKYGBgCgoKYGBge3J9Cm0xLjIgPC0gYnJtKGRheSB8IGNlbnMoY2VucykgfiAwICsgdGVtcHNjLAogICAgICAgICAgICBmYW1pbHkgPSBleHBvbmVudGlhbCgpLAogICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiwgY2xhc3M9ImIiKSwKICAgICAgICAgICAgZGF0YSA9IGdlcm0uc3RkaSwgKQpgYGAKCmBgYHtyfQpzdW1tYXJ5KG0xLjIpCmBgYAoKCmBgYHtyfQpwcmVkaWN0KG0xLjIpICU+JSBoZWFkCmBgYApsb2cobWVhbiB0aW1lIHRvIGdlcm1pbmF0aW9uKQoKYGBge3J9CnByZWRpY3QobTEuMiwgbmV3ZGF0YSA9IGV4cGFuZF9ncmlkKHRlbXBzYz11bmlxdWUoZ2VybS5zdGRpJHRlbXBzYyksIGNlbnM9Yygibm9uZSIsICJyaWdodCIpKSkKYGBgCgpgYGB7cn0KcHJlY2lzKG0xLjEsIGRlcHRoID0gMikKYGBgCmFnYWluLCB0aGVzZSBhcmUgbG9nKG1lYW4gdGltZSB0byBnZXJtaW5hdGlvbikKCmxhbWJkYSBzaG91bGQgYmUgMS9tdSwgc28KCmBgYHtyfQojcGxvdCgxOjI4LHBleHAoMToyOCwgMS9leHAoMi40OCkpLCB0eXBlPSJsIiwgY29sPSJyZWQiKQojbGluZXMoMToyOCxwZXhwKDE6MjgsIDEvZXhwKDUuNjApKSwgdHlwZT0ibCIsIGNvbD0iYmx1ZSIpCiNjcmFzaGVzCmBgYAoKCmBgYHtyfQpnZXRfcHJpb3IoIGRheSB8IGNlbnMoY2VucykgfiAwICsgdGVtcHNjLAogICAgICAgICAgIGZhbWlseSA9IEdhbW1hKCksCiAgICAgICAgICAgZGF0YSA9IGdlcm0uc3RkaSkKYGBgCgoKYGBge3J9Cm0xLjMgPC0gYnJtKGRheSB8IGNlbnMoY2VucykgfiAwICsgdGVtcHNjLAogICAgICAgICAgICBmYW1pbHkgPSBHYW1tYSgpLAogICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiwgY2xhc3M9ImIiLCBsYj0wKSwKICAgICAgICAgICAgZGF0YSA9IGdlcm0uc3RkaSkKYGBgCgpgYGB7cn0Kc3VtbWFyeShtMS4zKQpgYGAKCmBgYHtyfQpuZXdkYXRhIDwtIGV4cGFuZF9ncmlkKHRlbXBzYz11bmlxdWUoZ2VybS5zdGRpJHRlbXBzYyksIGNlbnM9Yygibm9uZSIsICJyaWdodCIpKQpjYmluZChuZXdkYXRhLCBwcmVkaWN0KG0xLjMsIG5ld2RhdGEgPSBuZXdkYXRhKSkKYGBgCgpUaGUgcHJvYmxlbSB3aXRoIHRoZXNlIG1vZGVscyBpcyB0aGF0IHRoZXkgYXNzdW1lIHRoYXQgdGhlIGNlbnNvcmVkIG9ic2VydmF0aW9ucyBhcmUgY29udHJvbGxlZCBieSB0aGUgc2FtZSBwcm9jZXNzIHRoYXQgaXMgY29udHJvbGxpbmcgdGhlIG9uZXMgdGhhdCBkbyBnZXJtaW5hdGUuICBJIGRvIG5vdCB0aGluayB0aGF0IGlzIGNvcnJlY3QsIHdlIG5lZWQgYSBtaXh0dXJlIG1vZGVsLgoKIyMgVHJ5aW5nIGEgbWl4dHVyZSBtb2RlbCBmb3IgYSBaSSBleHBvbmVudGlhbAoKbG9vayBhdCB1bGFtIGNlbnNvcmVkIGV4cG9uZW50aWFsCmBgYHtyfQpyZXRoaW5raW5nOjpzdGFuY29kZShtMS4xKQpgYGAKCmxvb2sgYXQgYnJtcyBjZW5zb3JlZCBleHBvbmVudGlhbApgYGB7cn0Kc3RhbmNvZGUobTEuMikKYGBgCgpaZXJvIGluZmxhdGVkLCBhcyBnZW5lcmF0ZWQgYnkgdWxhbToKYGBge3J9CiAKIyBkZWZpbmUgcGFyYW1ldGVycwpwcm9iX2RyaW5rIDwtIDAuMiAjIDIwJSBvZiBkYXlzCnJhdGVfd29yayA8LSAxICMgYXZlcmFnZSAxIG1hbnVzY3JpcHQgcGVyIGRheQojIHNhbXBsZSBvbmUgeWVhciBvZiBwcm9kdWN0aW9uIApOIDwtIDM2NQojIHNpbXVsYXRlIGRheXMgbW9ua3MgZHJpbmsgCnNldC5zZWVkKDM2NSkKZHJpbmsgPC0gcmJpbm9tKCBOICwgMSAsIHByb2JfZHJpbmsgKQojIHNpbXVsYXRlIG1hbnVzY3JpcHRzIGNvbXBsZXRlZAp5IDwtICgxLWRyaW5rKSpycG9pcyggTiAsIHJhdGVfd29yayApCgpyZXRoaW5raW5nOjpzdGFuY29kZSh1bGFtKCBhbGlzdCgKICB5fHk+MCB+IGN1c3RvbSggbG9nMW0ocCkgKyBwb2lzc29uX2xwbWYoeXxsYW1iZGEpICksCiAgeXx5PT0wIH4gY3VzdG9tKCBsb2dfbWl4KCBwICwgMCAsIHBvaXNzb25fbHBtZigwfGxhbWJkYSkgKSApLAogIGxvZ2l0KHApIDwtIGFwLAogIGxvZyhsYW1iZGEpIDwtIGFsLAogIGFwIH4gZG5vcm0oLTEuNSwxKSwKICBhbCB+IGRub3JtKDEsMC41KSksCiAgZGF0YT1saXN0KHk9YXMuaW50ZWdlcih5KSkgLCBzYW1wbGU9RkFMU0UpKQpgYGAKClplcm8gaW5mbGF0ZWQgZXhwb25lbnRpYWwgbW9kZWwsIGF0dGVtcHQgMQpgYGB7cn0KZCA8LSBsaXN0KE49bnJvdyhnZXJtLnN0ZGkpLAogICAgICAgICAgZ2VybT1nZXJtLnN0ZGkkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgpzdGFubW9kZWwxLjQgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhcyBmb3IgdGhlIGxpbmVhciBtb2RlbAogICAgcmVhbCBhcDsgLy8gcHJvYmFiaWxpdHkgb2YgZHJpbmsKfQptb2RlbHsKICAgIHJlYWwgcDsKICAgIHZlY3RvcltOXSBsYW1iZGE7CiAgICB2ZWN0b3JbTl0gbXU7CiAgICBhIH4gbm9ybWFsKCAwICwgMSApOwogICAgYXAgfiBub3JtYWwoIDAsIDEgKTsKICAgIHAgPSBhcDsKICAgIHAgPSBpbnZfbG9naXQocCk7IC8vIGludmVyc2UgbGluayBmdW5jdGlvbiBmb3IgcHJvYmFiaWxpdHkgZHJpbmsKICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07IC8vIGxpbmVhciBtb2RlbCBmb3IgbXUKICAgICAgICBtdVtpXSA9IGV4cChtdVtpXSk7ICAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24gcGFydCAxCiAgICAgICAgbGFtYmRhW2ldID0gMS9tdVtpXTsgLy8gaW52ZXJzZSBsaW5rIGZ1bmN0aW9uIHBhcnQgMgogICAgfQoKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDAgKSB0YXJnZXQgKz0gbG9nX21peChwLCAwLCBleHBvbmVudGlhbF9sY2NkZihkYXlbaV0gfCBsYW1iZGFbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSB0YXJnZXQgKz0gbG9nMW0ocCkgKyAgZXhwb25lbnRpYWxfbHBkZihkYXlbaV0gfCBsYW1iZGFbaV0pOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGV4cG9uZW50aWFsKCBsYW1iZGFbaV0gKTsKfQoiCm0xLjQgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuNCwgZGF0YT1kKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuNCwgZGVwdGggPSAyKQpgYGAKTW9kZWwgMS40IGZpdHMgYSBzaW5nbGUgZG9ybWFuY3kgcmF0ZSBmb3IgYWxsLCBidXQgcHJvYmFibHkgd2UgbmVlZCB0byBoYXZlIHNlcGFyYXRlIGRvcm1hbmN5IHJhdGUgZm9yIGVhY2ggdGVtcAoKdGFrZSBhIGxvb2sgYXQgc3RhbmNvZGUgZm9yIGEgYmlub21pYWwgbW9kZWwgdG8gc2VlIHdoYXQgdGhpcyBsb29rcyBsaWtlCgpgYGB7cn0KZGIgPC0gbGlzdChjdW11bGF0aXZlX2dlcm09Z2VybSRjdW11bGF0aXZlX2dlcm0sIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0kdGVtcHMpKSkKCnJldGhpbmtpbmc6OnN0YW5jb2RlKHVsYW0oYWxpc3QoCiAgY3VtdWxhdGl2ZV9nZXJtIH4gZGJlcm4obXUpLAogIG11IDwtIGxvZ2l0KGEgKyBiW3RlbXBzXSksCiAgYSB+IGRub3JtKDAsMSksCiAgYlt0ZW1wc10gfiBkbm9ybSgwLDEpKSwKICBkYXRhID0gZGIsIHNhbXBsZT1GQUxTRSwgbG9nX2xpayA9IFRSVUUpKQpgYGAKVGhlcmUgaXMgbm8gaW52ZXJzZSBsb2dpdCBiZWNhdXNlIGl0IGlzIGJ1aWx0IGluIHRvIHRoZSBiZXJub3VsbGkgZnVuY3Rpb24KClRyeSB0byB3cml0ZSBhIG1vZGVsIHdpdGggZGlmZmVyZW50IGRvcm1hbmN5IHByb2JzIGZvciB0aGUgZGlmZmVyZW50IHRlbXBzLgpgYGB7cn0KZCA8LSBsaXN0KE49bnJvdyhnZXJtLnN0ZGkpLAogICAgICAgICAgZ2VybT1nZXJtLnN0ZGkkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgpzdGFubW9kZWwxLjUgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgZXhwb25lbnRpYWwgY3VydmUsIG9uZSBmb3IgZWFjaCB0ZW1wCiAgICB2ZWN0b3JbOF0gYXA7IC8vIGFscGhhIGZvciB0aGUgcHJvcG9ydGlvbiBkb3JtYW50LCBvbmUgZm9yIGVhY2ggdGVtcAp9Cm1vZGVsewogICAgdmVjdG9yW05dIHA7CiAgICB2ZWN0b3JbTl0gbGFtYmRhOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIDEgKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAxICk7CiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICAgICAgbXVbaV0gPSBleHAobXVbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIGxhbWJkYVtpXSA9IDEvbXVbaV07CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZXhwb25lbnRpYWxfbGNjZGYoZGF5W2ldIHwgbGFtYmRhW2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgdGFyZ2V0ICs9IGxvZzFtKHBbaV0pICsgIGV4cG9uZW50aWFsX2xwZGYoZGF5W2ldIHwgbGFtYmRhW2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBleHBvbmVudGlhbCggbGFtYmRhW2ldICk7Cn0KIgptMS41IDwtIHN0YW4obW9kZWxfY29kZT1zdGFubW9kZWwxLjUsIGRhdGE9ZCkKYGBgCgpgYGB7cn0KcHJlY2lzKG0xLjUsIGRlcHRoID0gMikKYGBgCgpgYGB7cn0Kc2hpbnlzdGFuOjpsYXVuY2hfc2hpbnlzdGFuKG0xLjUpCmBgYAoKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuNSkKaGVhZChwb3N0KQpgYGAKCmBgYHtyfQpwb3N0LmNvciA8LSBjb3IocG9zdCkgJT4lIHJvdW5kKDIpCnBvc3QuY29yW2Ficyhwb3N0LmNvcik8LjI1XSA8LSBOQQpwb3N0LmNvcgpgYGAKCgpub3cgbmVlZCB0byBtYWtlIHByZWRpY3Rpb25zCgpgYGB7cn0KcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYVsiCikpICMgdGhlc2UgYXJlIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9sYW1iZGEgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBzdW1tYXJpc2VfYWxsKG1lYW4pICU+JSAjIHBvc3RlcmlvciBtZWFuCiAgbXV0YXRlX2FsbCh+IDEvZXhwKC4pKSAjIGNvbnZlcnQgdG8gbGFtYmRhIChpbnZlcnNlIGxpbmspCm11X2xhbWJkYQpgYGAKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCmBgYAoKYGBge3J9CnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZSh0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICBsYW1iZGE9dChtdV9sYW1iZGEpLAogICAgICAgICAgICAgICAgICAgICAgICBwPXQobXVfcCkpCmBgYAoKcGxvdCwgaW5jbHVkZSBvcmlnaW5hbCBkYXRhCmBgYHtyfQpwb3N0X3Bsb3QgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBleHAoZGF5LCByYXRlPWxhbWJkYSkgKiAoMS1wKSkKCnN0ZGkucGxvdCA8LSBnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiKSAlPiUgCiAgc2VsZWN0KGRheSwgdGVtcHMsIGN1bXVsYXRpdmVfZ2VybSwgdG90YWxfc2VlZHMpICU+JQogIG11dGF0ZSh0ZW1wcz1hcy5mYWN0b3IodGVtcHMpLAogICAgICAgICBwcm9wX2dlcm09Y3VtdWxhdGl2ZV9nZXJtL3RvdGFsX3NlZWRzKQoKcG9zdF9wbG90ICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGRhdGE9c3RkaS5wbG90KQogIApgYGAKCmFkanVzdCBwcmlvcnMKCmBgYHtyfQpkIDwtIGxpc3QoTj1ucm93KGdlcm0uc3RkaSksCiAgICAgICAgICBnZXJtPWdlcm0uc3RkaSRnZXJtLAogICAgICAgICAgdGVtcHM9YXMubnVtZXJpYyhhcy5mYWN0b3IoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICBkYXk9Z2VybS5zdGRpJGRheSkKCnN0YW5tb2RlbDEuNWEgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgZXhwb25lbnRpYWwgY3VydmUsIG9uZSBmb3IgZWFjaCB0ZW1wCiAgICB2ZWN0b3JbOF0gYXA7IC8vIGFscGhhIGZvciB0aGUgcHJvcG9ydGlvbiBkb3JtYW50LCBvbmUgZm9yIGVhY2ggdGVtcAp9Cm1vZGVsewogICAgdmVjdG9yW05dIHA7CiAgICB2ZWN0b3JbTl0gbGFtYmRhOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIDIgKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAyICk7CiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICAgICAgbXVbaV0gPSBleHAobXVbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIGxhbWJkYVtpXSA9IDEvbXVbaV07CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZXhwb25lbnRpYWxfbGNjZGYoZGF5W2ldIHwgbGFtYmRhW2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgdGFyZ2V0ICs9IGxvZzFtKHBbaV0pICsgIGV4cG9uZW50aWFsX2xwZGYoZGF5W2ldIHwgbGFtYmRhW2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBleHBvbmVudGlhbCggbGFtYmRhW2ldICk7Cn0KIgptMS41YSA8LSBzdGFuKG1vZGVsX2NvZGU9c3Rhbm1vZGVsMS41YSwgZGF0YT1kKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuNWEsIGRlcHRoID0gMikKYGBgCgpgYGB7cn0Kc2hpbnlzdGFuOjpsYXVuY2hfc2hpbnlzdGFuKG0xLjVhKQpgYGAKCgpgYGB7cn0KcG9zdCA8LSBhcy5kYXRhLmZyYW1lKG0xLjVhKQoKcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYVsiCikpICMgdGhlc2UgYXJlIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9sYW1iZGEgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBzdW1tYXJpc2VfYWxsKG1lYW4pICU+JSAjIHBvc3RlcmlvciBtZWFuCiAgbXV0YXRlX2FsbCh+IDEvZXhwKC4pKSAjIGNvbnZlcnQgdG8gbGFtYmRhIChpbnZlcnNlIGxpbmspCm11X2xhbWJkYQoKcG9zdF9sb2dpdHAgPC0gcG9zdCAlPiUgc2VsZWN0KHN0YXJ0c193aXRoKCJhcCIpKSAjbG9naXQgcCBkb3JtYW50Cm11X3AgPC0gcG9zdF9sb2dpdHAgJT4lCiAgc3VtbWFyaXplX2FsbChtZWFuKSAlPiUKICBtdXRhdGVfYWxsKGludl9sb2dpdCkKbXVfcAoKcG9zdGVyaW9yX2NvZWYgPC0gdGliYmxlKHRlbXBzPWFzLmZhY3Rvcih1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICAgICAgICAgICAgICAgIGxhbWJkYT10KG11X2xhbWJkYSksCiAgICAgICAgICAgICAgICAgICAgICAgIHA9dChtdV9wKSkKCnBvc3RfcGxvdCA8LSBleHBhbmRfZ3JpZChwb3N0ZXJpb3JfY29lZiwgZGF5PTE6MjgpICU+JQogIG11dGF0ZShwcm9wX2dlcm09cGV4cChkYXksIHJhdGU9bGFtYmRhKSAqICgxLXApKQoKc3RkaS5wbG90IDwtIGdlcm0gJT4lIGZpbHRlcihwb3BzPT0iU1RESSIpICU+JSAKICBzZWxlY3QoZGF5LCB0ZW1wcywgY3VtdWxhdGl2ZV9nZXJtLCB0b3RhbF9zZWVkcykgJT4lCiAgbXV0YXRlKHRlbXBzPWFzLmZhY3Rvcih0ZW1wcyksCiAgICAgICAgIHByb3BfZ2VybT1jdW11bGF0aXZlX2dlcm0vdG90YWxfc2VlZHMpCgpwb3N0X3Bsb3QgJT4lIAogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QpCiAgCmBgYApzdGlsbCBub3QgZ3JlYXQKCiMjIyBnYW1tYQoKTG9vayBhdCBnYW1tYSBzdGFuY29kZSBmcm9tIGJybXM6CgpgYGB7cn0KbWFrZV9zdGFuY29kZShkYXkgfiB0ZW1wc2MsIAogICAgICAgICAgICAgZmFtaWx5PSJHYW1tYSIsCiAgICAgICAgICAgICBwcmlvcj1zZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiksCiAgICAgICAgICAgICBkYXRhPWdlcm0uc3RkaSwKICAgICAgICAgICAgIHNhbXBsZT1GQUxTRSkKYGBgCgoKVHJ5IGRldmVsb3BpbmcgYSBaSSBnYW1tYSBtb2RlbDoKCnN0YW4gdXNlcyB0aGF0IGdhbW1hKGFscGhhLCBiZXRhKSBwYXJhbWV0ZXJpemF0aW9uIHdoZXJlICJhbHBoYSIgaXMgc2hhcGUiIGFuZCAiYmV0YSIgaXMgInJhdGUiLiAgU28gdHJ5IGhhdmluZyB0aGUgc2hhcGUgYmUgdGhlIHNhbWUgZm9yIGFsbCB0ZW1wcy4KYGBge3J9CmQgPC0gbGlzdChOPW5yb3coZ2VybS5zdGRpKSwKICAgICAgICAgIGdlcm09Z2VybS5zdGRpJGdlcm0sCiAgICAgICAgICB0ZW1wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMS43IDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHJlYWw8bG93ZXI9MD4gc2hhcGU7IC8vIHNob3VsZCBzZXQgbG93ZXIgYm91bmQKICAgIHZlY3Rvcls4XSBhOyAvLyBhbHBoYSBmb3IgdGhlIGdhbW1hIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIC41ICk7IC8vbmFycm93IHByaW9ycyB0byBvdmVyY29tZSBkaXZlcmdlbnQgdHJhbnNpdGlvbnMKICAgIGFwIH4gbm9ybWFsKCAwLCAxLjUgKTsgCiAgICBzaGFwZSB+IG5vcm1hbCgwLCAuNSk7IC8vbmFycm93IHByaW9ycyB0byBvdmVyY29tZSBkaXZlcmdlbnQgdHJhbnNpdGlvbnMKICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBbdGVtcHNbaV1dOwogICAgICAgIHBbaV0gPSBpbnZfbG9naXQocFtpXSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIHsKICAgICAgICBtdVtpXSA9IGFbdGVtcHNbaV1dOwogICAgICAgICAvLyBhcHBseSB0aGUgaW52ZXJzZSBsaW5rIGZ1bmN0aW9uCiAgICAgICAgbXVbaV0gPSBzaGFwZSAqIGV4cCgtKG11W2ldKSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDAgKSB0YXJnZXQgKz0gbG9nX21peChwW2ldLCAwLCBnYW1tYV9sY2NkZihkYXlbaV0gfCBzaGFwZSwgbXVbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArIGdhbW1hX2xwZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGdhbW1hKCBzaGFwZSwgbXVbaV0gKTsKfQoiCm0xLjcgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuNywgZGF0YT1kLCBjaGFpbnM9NCwgY29yZXM9NCwgY29udHJvbD1saXN0KGFkYXB0X2RlbHRhPS45OSkpCmBgYCAKCgpgYGB7cn0KcHJlY2lzKG0xLjcsIGRlcHRoID0gMikKYGBgCgoKYGBge3J9CnNoaW55c3Rhbjo6bGF1bmNoX3NoaW55c3RhbihtMS43KQpgYGAKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuNykKaGVhZChwb3N0KQpgYGAKCgpub3cgbmVlZCB0byBtYWtlIHByZWRpY3Rpb25zCgpgYGB7cn0KcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzaGFwZSwgc3RhcnRzX3dpdGgoImFbIgopKSAjIHRoZXNlIGFyZSB0aGUgc2hhcGUgYW5kIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9yYXRlIDwtIHBvc3RfbG9nYXZnZGF5cyAlPiUgCiAgbXV0YXRlKGFjcm9zcygtc2hhcGUsIH4gc2hhcGUqZXhwKC0oLikpKSApICU+JQogIHN1bW1hcml6ZShhY3Jvc3MoZXZlcnl0aGluZygpLCBtZWFuKSkgJT4lIHNlbGVjdCgtc2hhcGUpICMgcG9zdGVyaW9yIG1lYW4KbXVfcmF0ZQptdV9zaGFwZSA8LSBtZWFuKHBvc3Qkc2hhcGUpCm11X3NoYXBlCmBgYAoKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCmBgYAoKCnBsb3QgcHJlZGljdGlvbnMgYW5kIG9yaWdpbmFsIGRhdGEKYGBge3J9CnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZShzaGFwZT1tdV9zaGFwZSwgdGVtcHM9YXMuZmFjdG9yKHVuaXF1ZShnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgICAgICAgICAgICAgICAgcmF0ZT10KG11X3JhdGUpLAogICAgICAgICAgICAgICAgICAgICAgICBwPXQobXVfcCkpCgpwb3N0X3Bsb3QxLjcgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdDEuNyAlPiUgCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD10ZW1wcykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCkKICAKYGBgCgojIyMjIHBsYXlpbmcgd2l0aCBwcmlvcnM6CgpNYWtlIHNoYXBlIGV4cG9uZW50aWFsCmBgYHtyfQpkIDwtIGxpc3QoTj1ucm93KGdlcm0uc3RkaSksCiAgICAgICAgICBnZXJtPWdlcm0uc3RkaSRnZXJtLAogICAgICAgICAgdGVtcHM9YXMubnVtZXJpYyhhcy5mYWN0b3IoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICBkYXk9Z2VybS5zdGRpJGRheSkKCnN0YW5tb2RlbDEuN2EgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgcmVhbDxsb3dlcj0wPiBzaGFwZTsgLy8gc2hvdWxkIHNldCBsb3dlciBib3VuZAogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgZ2FtbWEgY3VydmUsIG9uZSBmb3IgZWFjaCB0ZW1wCiAgICB2ZWN0b3JbOF0gYXA7IC8vIGFscGhhIGZvciB0aGUgcHJvcG9ydGlvbiBkb3JtYW50LCBvbmUgZm9yIGVhY2ggdGVtcAp9Cm1vZGVsewogICAgdmVjdG9yW05dIHA7CiAgICB2ZWN0b3JbTl0gbXU7CiAgICBhIH4gbm9ybWFsKCAwICwgLjUgKTsgLy9uYXJyb3cgcHJpb3JzIHRvIG92ZXJjb21lIGRpdmVyZ2VudCB0cmFuc2l0aW9ucwogICAgYXAgfiBub3JtYWwoIDAsIDEuNSApOyAKICAgIHNoYXBlIH4gZXhwb25lbnRpYWwoLjUpOyAKICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBbdGVtcHNbaV1dOwogICAgICAgIHBbaV0gPSBpbnZfbG9naXQocFtpXSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIHsKICAgICAgICBtdVtpXSA9IGFbdGVtcHNbaV1dOwogICAgICAgICAvLyBhcHBseSB0aGUgaW52ZXJzZSBsaW5rIGZ1bmN0aW9uCiAgICAgICAgbXVbaV0gPSBzaGFwZSAqIGV4cCgtKG11W2ldKSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDAgKSB0YXJnZXQgKz0gbG9nX21peChwW2ldLCAwLCBnYW1tYV9sY2NkZihkYXlbaV0gfCBzaGFwZSwgbXVbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArIGdhbW1hX2xwZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGdhbW1hKCBzaGFwZSwgbXVbaV0gKTsKfQoiCm0xLjdhIDwtIHN0YW4obW9kZWxfY29kZT1zdGFubW9kZWwxLjdhLCBkYXRhPWQsIGNoYWlucz00LCBjb3Jlcz00LCBjb250cm9sPWxpc3QoYWRhcHRfZGVsdGE9Ljk5KSkKYGBgIAoKCmBgYHtyfQpwcmVjaXMobTEuN2EsIGRlcHRoID0gMikKYGBgCgoKYGBge3J9CnBvc3QgPC0gYXMuZGF0YS5mcmFtZShtMS43YSkKCnBvc3RfbG9nYXZnZGF5cyA8LSBwb3N0ICU+JSBzZWxlY3Qoc2hhcGUsIHN0YXJ0c193aXRoKCJhWyIKKSkgIyB0aGVzZSBhcmUgdGhlIHNoYXBlIGFuZCB0aGUgbG9nKG1lYW4oYXZnIHRpbWUgdG8gZ2VybWluYXRlKSkKbXVfcmF0ZSA8LSBwb3N0X2xvZ2F2Z2RheXMgJT4lIAogIG11dGF0ZShhY3Jvc3MoLXNoYXBlLCB+IHNoYXBlKmV4cCgtKC4pKSkgKSAlPiUKICBzdW1tYXJpemUoYWNyb3NzKGV2ZXJ5dGhpbmcoKSwgbWVhbikpICU+JSBzZWxlY3QoLXNoYXBlKSAjIHBvc3RlcmlvciBtZWFuCm11X3JhdGUKbXVfc2hhcGUgPC0gbWVhbihwb3N0JHNoYXBlKQoKcG9zdF9sb2dpdHAgPC0gcG9zdCAlPiUgc2VsZWN0KHN0YXJ0c193aXRoKCJhcCIpKSAjbG9naXQgcCBkb3JtYW50Cm11X3AgPC0gcG9zdF9sb2dpdHAgJT4lCiAgc3VtbWFyaXplX2FsbChtZWFuKSAlPiUKICBtdXRhdGVfYWxsKGludl9sb2dpdCkKbXVfcAoKcG9zdGVyaW9yX2NvZWYgPC0gdGliYmxlKHNoYXBlPW11X3NoYXBlLCB0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICByYXRlPXQobXVfcmF0ZSksCiAgICAgICAgICAgICAgICAgICAgICAgIHA9dChtdV9wKSkKCnBvc3RfcGxvdDEuN2EgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdDEuN2EgJT4lIAogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QpCiAgCmBgYApTaGFwZSBwYXJhbWV0ZXIgaXMgbGFyZ2VyIHRoYW4gaW4gMS43LiAgVGhlIGZpdCBpcyBzbGlnaHRseSBiZXR0ZXIgYXQgbGVhc3QgZm9yIDIwIGFuZCAyNS4KCm1ha2UgYSBjb21wYXJpc29uIHBsb3QKCmBgYHtyfQpwb3N0X3Bsb3RfY29tcGFyZSA8LSBiaW5kX3Jvd3Mob2xkPXBvc3RfcGxvdDEuNywgbmV3PXBvc3RfcGxvdDEuN2EsIC5pZCA9ICJmaXQiKSAKcG9zdF9wbG90X2NvbXBhcmUgJT4lCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD1zdHJfYyh0ZW1wcyxmaXQpLGxpbmV0eXBlPWZpdCkpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCwgYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzKSwgaW5oZXJpdC5hZXMgPSBGQUxTRSkKICAKYGBgCkp1c3QgYSBoYWlyIGJldHRlci4uLgoKY29udGludWUgdG8gd29yayB3aXRoIHByaW9ycwpgYGB7cn0KZCA8LSBsaXN0KE49bnJvdyhnZXJtLnN0ZGkpLAogICAgICAgICAgZ2VybT1nZXJtLnN0ZGkkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgZGF5PWdlcm0uc3RkaSRkYXkpCgpzdGFubW9kZWwxLjdiIDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHJlYWw8bG93ZXI9MD4gc2hhcGU7IC8vIHNob3VsZCBzZXQgbG93ZXIgYm91bmQKICAgIHZlY3Rvcls4XSBhOyAvLyBhbHBoYSBmb3IgdGhlIGdhbW1hIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIC43NSApOyAvLyBkaXZlcmdlbnQgd2l0aCAxIGFuZCB3aXRoIC43NQogICAgYXAgfiBub3JtYWwoIDAsIDEuNSApOyAKICAgIHNoYXBlIH4gZXhwb25lbnRpYWwoMS41KTsgLy8gZGl2ZXJnZW50IGlmIHRoaXMgaXMgLjUgYW5kIGEgaXMgMCwgMQogICAgZm9yIChpIGluIDE6TikgewogICAgICAgIHBbaV0gPSBhcFt0ZW1wc1tpXV07CiAgICAgICAgcFtpXSA9IGludl9sb2dpdChwW2ldKTsKICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICAgICAgIC8vIGFwcGx5IHRoZSBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgICAgICBtdVtpXSA9IHNoYXBlICogZXhwKC0obXVbaV0pKTsKICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMCApIHRhcmdldCArPSBsb2dfbWl4KHBbaV0sIDAsIGdhbW1hX2xjY2RmKGRheVtpXSB8IHNoYXBlLCBtdVtpXSkpOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgdGFyZ2V0ICs9IGxvZzFtKHBbaV0pICsgZ2FtbWFfbHBkZihkYXlbaV0gfCBzaGFwZSwgbXVbaV0pOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgZGF5W2ldIH4gZ2FtbWEoIHNoYXBlLCBtdVtpXSApOwp9CiIKbTEuN2IgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuN2IsIGRhdGE9ZCwgY2hhaW5zPTQsIGNvcmVzPTQsIGNvbnRyb2w9bGlzdChhZGFwdF9kZWx0YT0uOTkpKQpgYGAgCgoKYGBge3J9CnByZWNpcyhtMS43YiwgZGVwdGggPSAyKQpgYGAKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuN2IpCgpwb3N0X2xvZ2F2Z2RheXMgPC0gcG9zdCAlPiUgc2VsZWN0KHNoYXBlLCBzdGFydHNfd2l0aCgiYVsiCikpICMgdGhlc2UgYXJlIHRoZSBzaGFwZSBhbmQgdGhlIGxvZyhtZWFuKGF2ZyB0aW1lIHRvIGdlcm1pbmF0ZSkpCm11X3JhdGUgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBtdXRhdGUoYWNyb3NzKC1zaGFwZSwgfiBzaGFwZSpleHAoLSguKSkpICkgJT4lCiAgc3VtbWFyaXplKGFjcm9zcyhldmVyeXRoaW5nKCksIG1lYW4pKSAlPiUgc2VsZWN0KC1zaGFwZSkgIyBwb3N0ZXJpb3IgbWVhbgptdV9yYXRlCm11X3NoYXBlIDwtIG1lYW4ocG9zdCRzaGFwZSkKCnBvc3RfbG9naXRwIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYXAiKSkgI2xvZ2l0IHAgZG9ybWFudAptdV9wIDwtIHBvc3RfbG9naXRwICU+JQogIHN1bW1hcml6ZV9hbGwobWVhbikgJT4lCiAgbXV0YXRlX2FsbChpbnZfbG9naXQpCm11X3AKCnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZShzaGFwZT1tdV9zaGFwZSwgdGVtcHM9YXMuZmFjdG9yKHVuaXF1ZShnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgICAgICAgICAgICAgICAgcmF0ZT10KG11X3JhdGUpLAogICAgICAgICAgICAgICAgICAgICAgICBwPXQobXVfcCkpCgpwb3N0X3Bsb3QgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdCAlPiUgCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD10ZW1wcykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCkKICAKYGBgCm5vIGltcHJvdmVtZW50IG92ZXIgMS43YQoKIyMjIGJybXMgaHVyZGxlIGdhbW1hCgpgYGB7cn0KZ2V0X3ByaW9yKGJmKGRheSB+IDAgKyB0ZW1wc2MsIGh1IH4gMCArIHRlbXBzYyksCiAgICAgICAgICBmYW1pbHkgPSBodXJkbGVfZ2FtbWEsCiAgICAgICAgICBkYXRhPWdlcm0uc3RkaSkKYGBgCgpgYGB7cn0KbWFrZV9zdGFuY29kZShiZihkYXkgfiAwICsgdGVtcHNjLCBodSB+IDAgKyB0ZW1wc2MpLAogICAgICAgICAgZmFtaWx5ID0gaHVyZGxlX2dhbW1hLAogICAgICAgICAgcHJpb3I9c2V0X3ByaW9yKCJub3JtYWwoMCwxKSIpLAogICAgICAgICAgZGF0YT1nZXJtLnN0ZGkpCmBgYAoKCmBgYHtyfQptMS44IDwtIGJybShiZihkYXkgfiAwICsgdGVtcHNjLCBodSB+IDAgKyB0ZW1wc2MpLAogICAgICAgICAgZmFtaWx5ID0gaHVyZGxlX2dhbW1hLAogICAgICAgICAgcHJpb3I9YyhzZXRfcHJpb3IoIm5vcm1hbCgwLDEpIiksCiAgICAgICAgICAgICAgICAgIHNldF9wcmlvcigibm9ybWFsKDAsMSkiLCBkcGFyPSJodSIpKSwgIyBtdXN0IHNwZWNpZnkgdGhpcyEKICAgICAgICAgIGRhdGE9Z2VybS5zdGRpLCBjb3JlcyA9IDQpCmBgYAoKCmBgYHtyfQptMS44CmBgYAoKCmBgYHtyfQpwb3N0IDwtIGFzLmRhdGEuZnJhbWUobTEuOCkKaGVhZChwb3N0KQpgYGAKCndlaXJkIGZpdCEKCm5vdyBuZWVkIHRvIG1ha2UgcHJlZGljdGlvbnMKCmBgYHtyfQpwb3N0X2xvZ2F2Z2RheXMgPC0gcG9zdCAlPiUgc2VsZWN0KHNoYXBlLCBzdGFydHNfd2l0aCgiYl90ZW1wIgopKSAjIHRoZXNlIGFyZSB0aGUgc2hhcGUgYW5kIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9yYXRlIDwtIHBvc3RfbG9nYXZnZGF5cyAlPiUgCiAgbXV0YXRlKGFjcm9zcygtc2hhcGUsIH4gc2hhcGUqZXhwKC0oLikpKSApICU+JQogIHN1bW1hcml6ZShhY3Jvc3MoZXZlcnl0aGluZygpLCBtZWFuKSkgJT4lIHNlbGVjdCgtc2hhcGUpICMgcG9zdGVyaW9yIG1lYW4KbXVfcmF0ZQptdV9zaGFwZSA8LSBtZWFuKHBvc3Qkc2hhcGUpCmBgYAoKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImJfaHUiKSkgI2xvZ2l0IHAgZG9ybWFudAptdV9wIDwtIHBvc3RfbG9naXRwICU+JQogIHN1bW1hcml6ZV9hbGwobWVhbikgJT4lCiAgbXV0YXRlX2FsbChpbnZfbG9naXQpCm11X3AKYGBgCgpgYGB7cn0KcG9zdGVyaW9yX2NvZWYgPC0gdGliYmxlKHNoYXBlPW11X3NoYXBlLCB0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICByYXRlPXQobXVfcmF0ZSksCiAgICAgICAgICAgICAgICAgICAgICAgIHA9dChtdV9wKSkKYGBgCgpwbG90IHByZWRpY3Rpb25zIGFuZCBvcmlnaW5hbCBkYXRhCmBgYHtyfQoKCnN0ZGkucGxvdCA8LSBnZXJtICU+JSBmaWx0ZXIocG9wcz09IlNUREkiKSAlPiUgCiAgc2VsZWN0KGRheSwgdGVtcHMsIGN1bXVsYXRpdmVfZ2VybSwgdG90YWxfc2VlZHMpICU+JQogIG11dGF0ZSh0ZW1wcz1hcy5mYWN0b3IodGVtcHMpLAogICAgICAgICBwcm9wX2dlcm09Y3VtdWxhdGl2ZV9nZXJtL3RvdGFsX3NlZWRzKQoKcG9zdF9wbG90ICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGRhdGE9c3RkaS5wbG90KQogIApgYGAKCiMjIGFsbG93IGJvdGggc2hhcGUgYW5kIHJhdGUgb2YgWkkgZ2FtbWEgdG8gYmUgcHJlZGljdGVkIGJ5IHRlbXAKCkRvbid0IHJ1bi4gIERvZXNuJ3Qgc2FtcGxlIHdlbGwgYmVjYXVzZSBzaGFwZSBhbmQgcmF0ZSBhcmUgY29ycmVsYXRlZC4KYGBge3IsIGV2YWw9RkFMU0V9CmQgPC0gbGlzdChOPW5yb3coZ2VybS5zdGRpKSwKICAgICAgICAgIGdlcm09Z2VybS5zdGRpJGdlcm0sCiAgICAgICAgICB0ZW1wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMS45IDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHZlY3Rvcls4XSBhczsgLy8gYWxwaGEgZm9yIGdhbWEgc2hhcGUKICAgIHZlY3Rvcls4XSBhcjsgLy8gYWxwaGEgZm9yIHRoZSBnYW1tYSByYXRlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIHNoYXBlOwogICAgdmVjdG9yW05dIHJhdGU7CiAgICBhcyB+IG5vcm1hbCggMCwgMSApOwogICAgYXIgfiBub3JtYWwoIDAsIDEgKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAxLjUgKTsgCiAgICAKICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBbdGVtcHNbaV1dOwogICAgICAgIHBbaV0gPSBpbnZfbG9naXQocFtpXSk7CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIHsKICAgICAgICBzaGFwZVtpXSA9IGFzW3RlbXBzW2ldXTsKICAgICAgICByYXRlW2ldID0gYXJbdGVtcHNbaV1dOwogICAgICAgICAvLyBhcHBseSB0aGUgaW52ZXJzZSBsaW5rIGZ1bmN0aW9ucwogICAgICAgIHNoYXBlW2ldID0gZXhwKHNoYXBlW2ldKTsKICAgICAgICByYXRlW2ldID0gZXhwKHJhdGVbaV0pOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZ2FtbWFfbGNjZGYoZGF5W2ldIHwgc2hhcGVbaV0sIHJhdGVbaV0pKTsKICAgIGZvciAoIGkgaW4gMTpOICkgCiAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArIGdhbW1hX2xwZGYoZGF5W2ldIHwgc2hhcGVbaV0sIHJhdGVbaV0pOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAxICkgZGF5W2ldIH4gZ2FtbWEoIHNoYXBlW2ldLCByYXRlIFtpXSApOwp9CmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKCn0KIgptMS45IDwtIHN0YW4obW9kZWxfY29kZT1zdGFubW9kZWwxLjksIGRhdGE9ZCwgY2hhaW5zPTQsIGNvcmVzPTQsIGNvbnRyb2w9bGlzdChhZGFwdF9kZWx0YT0uOTkpKQpgYGAgCgoKYGBge3J9CnByZWNpcyhtMS45LCBkZXB0aCA9IDIpCmBgYAoKYGBge3J9CnBvc3QgPC0gYXMuZGF0YS5mcmFtZShtMS45KQoKcG9zdF9sb2dhdmdkYXlzIDwtIHBvc3QgJT4lIHNlbGVjdChzaGFwZSwgc3RhcnRzX3dpdGgoImFbIgopKSAjIHRoZXNlIGFyZSB0aGUgc2hhcGUgYW5kIHRoZSBsb2cobWVhbihhdmcgdGltZSB0byBnZXJtaW5hdGUpKQptdV9yYXRlIDwtIHBvc3RfbG9nYXZnZGF5cyAlPiUgCiAgbXV0YXRlKGFjcm9zcygtc2hhcGUsIH4gc2hhcGUqZXhwKC0oLikpKSApICU+JQogIHN1bW1hcml6ZShhY3Jvc3MoZXZlcnl0aGluZygpLCBtZWFuKSkgJT4lIHNlbGVjdCgtc2hhcGUpICMgcG9zdGVyaW9yIG1lYW4KbXVfcmF0ZQptdV9zaGFwZSA8LSBtZWFuKHBvc3Qkc2hhcGUpCgpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCgpwb3N0ZXJpb3JfY29lZiA8LSB0aWJibGUoc2hhcGU9bXVfc2hhcGUsIHRlbXBzPWFzLmZhY3Rvcih1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICAgICAgICAgICAgICAgIHJhdGU9dChtdV9yYXRlKSwKICAgICAgICAgICAgICAgICAgICAgICAgcD10KG11X3ApKQoKcG9zdF9wbG90IDwtIGV4cGFuZF9ncmlkKHBvc3Rlcmlvcl9jb2VmLCBkYXk9MToyOCkgJT4lCiAgbXV0YXRlKHByb3BfZ2VybT1wZ2FtbWEoZGF5LCBzaGFwZT1tdV9zaGFwZSwgcmF0ZT1yYXRlKSAqICgxLXApKQoKc3RkaS5wbG90IDwtIGdlcm0gJT4lIGZpbHRlcihwb3BzPT0iU1RESSIpICU+JSAKICBzZWxlY3QoZGF5LCB0ZW1wcywgY3VtdWxhdGl2ZV9nZXJtLCB0b3RhbF9zZWVkcykgJT4lCiAgbXV0YXRlKHRlbXBzPWFzLmZhY3Rvcih0ZW1wcyksCiAgICAgICAgIHByb3BfZ2VybT1jdW11bGF0aXZlX2dlcm0vdG90YWxfc2VlZHMpCgpwb3N0X3Bsb3QgJT4lIAogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QpCiAgCmBgYAoKIyMgbG9nIG5vcm1hbAoKYGBge3J9Cm1ha2Vfc3RhbmNvZGUoZGF5IH4gdGVtcHNjLCBmYW1pbHkgPSBsb2dub3JtYWwoKSwgZGF0YT1nZXJtLnN0ZGkpCmBgYAoKCmBgYHtyfQpkIDwtIGxpc3QoTj1ucm93KGdlcm0uc3RkaSksCiAgICAgICAgICBnZXJtPWdlcm0uc3RkaSRnZXJtLAogICAgICAgICAgdGVtcHM9YXMubnVtZXJpYyhhcy5mYWN0b3IoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICBkYXk9Z2VybS5zdGRpJGRheSkKCnN0YW5tb2RlbDEuMTAgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgdmVjdG9yWzhdIGE7IC8vIGFscGhhIGZvciB0aGUgbG9nbm9ybWFsIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKICAgIHJlYWw8bG93ZXI9MD4gc2lnbWE7Cn0KbW9kZWx7CiAgICB2ZWN0b3JbTl0gcDsKICAgIHZlY3RvcltOXSBtdTsKICAgIGEgfiBub3JtYWwoIDAgLCAyICk7CiAgICBzaWdtYSB+IGV4cG9uZW50aWFsICgxKTsKICAgIGFwIH4gbm9ybWFsKCAwLCAyICk7CiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOyAvLyBpbnZlcnNlIGxpbmsgZnVuY3Rpb24KICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYVt0ZW1wc1tpXV07CiAgICB9CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgbG9nbm9ybWFsX2xjY2RmKGRheVtpXSB8IG11W2ldLCBzaWdtYSkpOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIHRhcmdldCArPSBsb2cxbShwW2ldKSArICBsb2dub3JtYWxfbHBkZihkYXlbaV0gfCBtdVtpXSwgc2lnbWEpOwogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgICBpZiAoIGdlcm1baV0gPT0gMSApIGRheVtpXSB+IGxvZ25vcm1hbCggbXVbaV0sIHNpZ21hICk7Cn0KIgptMS4xMCA8LSBzdGFuKG1vZGVsX2NvZGU9c3Rhbm1vZGVsMS4xMCwgZGF0YT1kKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuMTAsIGRlcHRoID0gMikKYGBgCgpgYGB7cn0Kc2hpbnlzdGFuOjpsYXVuY2hfc2hpbnlzdGFuKG0xLjEwKQpgYGAKCgpgYGB7cn0KcG9zdCA8LSBhcy5kYXRhLmZyYW1lKG0xLjEwKQpoZWFkKHBvc3QpCmBgYAoKbm93IG5lZWQgdG8gbWFrZSBwcmVkaWN0aW9ucwoKYGBge3J9CnBvc3RfbG9nYXZnZGF5cyA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFbIgopKSAjIHRoZXNlIGFyZSB0aGUgbG9nKG1lYW4oYXZnIHRpbWUgdG8gZ2VybWluYXRlKSkKbXUgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBzdW1tYXJpc2VfYWxsKG1lYW4pIAoKbXVfc2lnbWEgPC0gbWVhbihwb3N0JHNpZ21hKQpgYGAKCmBgYHtyfQpwb3N0X2xvZ2l0cCA8LSBwb3N0ICU+JSBzZWxlY3Qoc3RhcnRzX3dpdGgoImFwIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCmBgYAoKYGBge3J9CnBvc3Rlcmlvcl9jb2VmIDwtIHRpYmJsZSh0ZW1wcz1hcy5mYWN0b3IodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgICAgICAgICAgICAgICBtdT10KG11KSwKICAgICAgICAgICAgICAgICAgICAgICAgcD10KG11X3ApLAogICAgICAgICAgICAgICAgICAgICAgICBzaWdtYT1tdV9zaWdtYSkKYGBgCgpwbG90LCBpbmNsdWRlIG9yaWdpbmFsIGRhdGEKYGBge3J9CnBvc3RfcGxvdCA8LSBleHBhbmRfZ3JpZChwb3N0ZXJpb3JfY29lZiwgZGF5PTE6MjgpICU+JQogIG11dGF0ZShwcm9wX2dlcm09cGxub3JtKGRheSwgbWVhbmxvZz1tdSwgc2Rsb2c9c2lnbWEpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdCAlPiUgCiAgZ2dwbG90KGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyxncm91cD10ZW1wcykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChkYXRhPXN0ZGkucGxvdCkKICAKYGBgCgojIHdlZWsyCgojIyBtdiBub3JtIGZvciByYXRlIGFuZCBwPwoKTG9vayBhdCBTdGFuY29kZSBmcm9tIGNoYXB0ZXIgMTQgbXZub3JtCmBgYHtyfQojIyBSIGNvZGUgMTQuMQphIDwtIDMuNSAgICAgICAgICAgICMgYXZlcmFnZSBtb3JuaW5nIHdhaXQgdGltZQpiIDwtICgtMSkgICAgICAgICAgICMgYXZlcmFnZSBkaWZmZXJlbmNlIGFmdGVybm9vbiB3YWl0IHRpbWUKc2lnbWFfYSA8LSAxICAgICAgICAjIHN0ZCBkZXYgaW4gaW50ZXJjZXB0cwpzaWdtYV9iIDwtIDAuNSAgICAgICMgc3RkIGRldiBpbiBzbG9wZXMKcmhvIDwtICgtMC43KSAgICAgICAjIGNvcnJlbGF0aW9uIGJldHdlZW4gaW50ZXJjZXB0cyBhbmQgc2xvcGVzCgojIyBSIGNvZGUgMTQuMgpNdSA8LSBjKCBhICwgYiApCgojIyBSIGNvZGUgMTQuMwpjb3ZfYWIgPC0gc2lnbWFfYSpzaWdtYV9iKnJobwpTaWdtYSA8LSBtYXRyaXgoIGMoc2lnbWFfYV4yLGNvdl9hYixjb3ZfYWIsc2lnbWFfYl4yKSAsIG5jb2w9MiApCgoKIyMgUiBjb2RlIDE0LjUKc2lnbWFzIDwtIGMoc2lnbWFfYSxzaWdtYV9iKSAjIHN0YW5kYXJkIGRldmlhdGlvbnMKUmhvIDwtIG1hdHJpeCggYygxLHJobyxyaG8sMSkgLCBucm93PTIgKSAjIGNvcnJlbGF0aW9uIG1hdHJpeAoKIyBub3cgbWF0cml4IG11bHRpcGx5IHRvIGdldCBjb3ZhcmlhbmNlIG1hdHJpeApTaWdtYSA8LSBkaWFnKHNpZ21hcykgJSolIFJobyAlKiUgZGlhZyhzaWdtYXMpCgojIyBSIGNvZGUgMTQuNgpOX2NhZmVzIDwtIDIwCgojIyBSIGNvZGUgMTQuNwpzZXQuc2VlZCg1KSAjIHVzZWQgdG8gcmVwbGljYXRlIGV4YW1wbGUKdmFyeV9lZmZlY3RzIDwtIE1BU1M6Om12cm5vcm0oIE5fY2FmZXMgLCBNdSAsIFNpZ21hICkKCiMjIFIgY29kZSAxNC44CmFfY2FmZSA8LSB2YXJ5X2VmZmVjdHNbLDFdCmJfY2FmZSA8LSB2YXJ5X2VmZmVjdHNbLDJdCgojIyBSIGNvZGUgMTQuMTAKc2V0LnNlZWQoMjIpCk5fdmlzaXRzIDwtIDEwCmFmdGVybm9vbiA8LSByZXAoMDoxLE5fdmlzaXRzKk5fY2FmZXMvMikKY2FmZV9pZCA8LSByZXAoIDE6Tl9jYWZlcyAsIGVhY2g9Tl92aXNpdHMgKQptdSA8LSBhX2NhZmVbY2FmZV9pZF0gKyBiX2NhZmVbY2FmZV9pZF0qYWZ0ZXJub29uCnNpZ21hIDwtIDAuNSAgIyBzdGQgZGV2IHdpdGhpbiBjYWZlcwp3YWl0IDwtIHJub3JtKCBOX3Zpc2l0cypOX2NhZmVzICwgbXUgLCBzaWdtYSApCmQgPC0gZGF0YS5mcmFtZSggY2FmZT1jYWZlX2lkICwgYWZ0ZXJub29uPWFmdGVybm9vbiAsIHdhaXQ9d2FpdCApCgoKcmV0aGlua2luZzo6c3RhbmNvZGUodWxhbSgKICAgIGFsaXN0KAogICAgICAgIHdhaXQgfiBub3JtYWwoIG11ICwgc2lnbWEgKSwKICAgICAgICBtdSA8LSBhX2NhZmVbY2FmZV0gKyBiX2NhZmVbY2FmZV0qYWZ0ZXJub29uLAogICAgICAgIGMoYV9jYWZlLGJfY2FmZSlbY2FmZV0gfiBtdWx0aV9ub3JtYWwoIGMoYSxiKSAsIFJobyAsIHNpZ21hX2NhZmUgKSwKICAgICAgICBhIH4gbm9ybWFsKDUsMiksCiAgICAgICAgYiB+IG5vcm1hbCgtMSwwLjUpLAogICAgICAgIHNpZ21hX2NhZmUgfiBleHBvbmVudGlhbCgxKSwKICAgICAgICBzaWdtYSB+IGV4cG9uZW50aWFsKDEpLAogICAgICAgIFJobyB+IGxral9jb3JyKDIpCiAgICApICwgZGF0YT1kICwgc2FtcGxlID0gRkFMU0UgKSkKYGBgCgoKCk1ha2Ugc2hhcGUgZXhwb25lbnRpYWwKYGBge3J9CmQgPC0gbGlzdChOPW5yb3coZ2VybS5zdGRpKSwKICAgICAgICAgIGdlcm09Z2VybS5zdGRpJGdlcm0sCiAgICAgICAgICB0ZW1wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtLnN0ZGkkdGVtcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMS43bXYgPC0KICAiCmRhdGF7CiAgICBpbnQ8bG93ZXI9MT4gTjsgIC8vIG51bWJlciBvZiBvYnNlcnZhdGlvbnMKICAgIGludCBnZXJtW05dOwogICAgdmVjdG9yW05dIGRheTsKICAgIGludCB0ZW1wc1tOXTsKfQpwYXJhbWV0ZXJzewogICAgcmVhbDxsb3dlcj0wPiBzaGFwZTsgLy8gc2hvdWxkIHNldCBsb3dlciBib3VuZAogICAgdmVjdG9yWzhdIGFfdGVtcDsgLy8gYWxwaGEgZm9yIHRoZSBnYW1tYSBjdXJ2ZSwgb25lIGZvciBlYWNoIHRlbXAKICAgIHZlY3Rvcls4XSBhcF90ZW1wOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKICAgIHJlYWwgYTsKICAgIHJlYWwgYXA7CiAgICB2ZWN0b3I8bG93ZXI9MD5bMl0gc2lnbWFfdGVtcDsKICAgIGNvcnJfbWF0cml4WzJdIFJobzsKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgUmhvIH4gbGtqX2NvcnIoMik7CiAgICBzaWdtYV90ZW1wIH4gZXhwb25lbnRpYWwoMSk7CiAgICAKICAgIGEgfiBub3JtYWwoIDAgLCAuNSApOyAvL25hcnJvdyBwcmlvcnMgdG8gb3ZlcmNvbWUgZGl2ZXJnZW50IHRyYW5zaXRpb25zCiAgICBhcCB+IG5vcm1hbCggMCwgMS41ICk7IAogICAgc2hhcGUgfiBleHBvbmVudGlhbCguNSk7IAp7CiAgICB2ZWN0b3JbMl0gWVlbOF07CiAgICB2ZWN0b3JbMl0gTVU7CiAgICBNVSA9IFthLCBhcF0nOwogICAgZm9yIChqIGluIDE6OCkgWVlbal0gPSBbYV90ZW1wW2pdLCBhcF90ZW1wW2pdXSc7CiAgICBZWSB+IG11bHRpX25vcm1hbChNVSwgcXVhZF9mb3JtX2RpYWcoUmhvLCBzaWdtYV90ZW1wKSk7Cn0KICAgIGZvciAoaSBpbiAxOk4pIHsKICAgICAgICBwW2ldID0gYXBfdGVtcFt0ZW1wc1tpXV07CiAgICAgICAgcFtpXSA9IGludl9sb2dpdChwW2ldKTsKICAgIH0KICAgIGZvciAoIGkgaW4gMTpOICkgewogICAgICAgIG11W2ldID0gYV90ZW1wW3RlbXBzW2ldXTsKICAgICAgICAgLy8gYXBwbHkgdGhlIGludmVyc2UgbGluayBmdW5jdGlvbgogICAgICAgIG11W2ldID0gc2hhcGUgKiBleHAoLShtdVtpXSkpOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZ2FtbWFfbGNjZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSB0YXJnZXQgKz0gbG9nMW0ocFtpXSkgKyBnYW1tYV9scGRmKGRheVtpXSB8IHNoYXBlLCBtdVtpXSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBnYW1tYSggc2hhcGUsIG11W2ldICk7Cn0KIgptMS43bXYgPC0gc3Rhbihtb2RlbF9jb2RlPXN0YW5tb2RlbDEuN212LCBkYXRhPWQsIGNoYWlucz00LCBjb3Jlcz00LCBjb250cm9sPWxpc3QoYWRhcHRfZGVsdGE9Ljk5KSkKYGBgIAoKMTAwIGRpdmVyZ2VudCB0cmFuc2l0aW9ucy4uLgoKYGBge3J9CnByZWNpcyhtMS43bXYsIGRlcHRoID0gMykKYGBgCgoKYGBge3J9CnBvc3QgPC0gYXMuZGF0YS5mcmFtZShtMS43bXYpCgpwb3N0X2xvZ2F2Z2RheXMgPC0gcG9zdCAlPiUgc2VsZWN0KHNoYXBlLCBzdGFydHNfd2l0aCgiYV8iCikpICMgdGhlc2UgYXJlIHRoZSBzaGFwZSBhbmQgdGhlIGxvZyhtZWFuKGF2ZyB0aW1lIHRvIGdlcm1pbmF0ZSkpCm11X3JhdGUgPC0gcG9zdF9sb2dhdmdkYXlzICU+JSAKICBtdXRhdGUoYWNyb3NzKC1zaGFwZSwgfiBzaGFwZSpleHAoLSguKSkpICkgJT4lCiAgc3VtbWFyaXplKGFjcm9zcyhldmVyeXRoaW5nKCksIG1lYW4pKSAlPiUgc2VsZWN0KC1zaGFwZSkgIyBwb3N0ZXJpb3IgbWVhbgptdV9yYXRlCm11X3NoYXBlIDwtIG1lYW4ocG9zdCRzaGFwZSkKCnBvc3RfbG9naXRwIDwtIHBvc3QgJT4lIHNlbGVjdChzdGFydHNfd2l0aCgiYXBfIikpICNsb2dpdCBwIGRvcm1hbnQKbXVfcCA8LSBwb3N0X2xvZ2l0cCAlPiUKICBzdW1tYXJpemVfYWxsKG1lYW4pICU+JQogIG11dGF0ZV9hbGwoaW52X2xvZ2l0KQptdV9wCgpwb3N0ZXJpb3JfY29lZiA8LSB0aWJibGUoc2hhcGU9bXVfc2hhcGUsIHRlbXBzPWFzLmZhY3Rvcih1bmlxdWUoZ2VybS5zdGRpJHRlbXBzKSksCiAgICAgICAgICAgICAgICAgICAgICAgIHJhdGU9dChtdV9yYXRlKSwKICAgICAgICAgICAgICAgICAgICAgICAgcD10KG11X3ApKQoKcG9zdF9wbG90MS43bXYgPC0gZXhwYW5kX2dyaWQocG9zdGVyaW9yX2NvZWYsIGRheT0xOjI4KSAlPiUKICBtdXRhdGUocHJvcF9nZXJtPXBnYW1tYShkYXksIHNoYXBlPW11X3NoYXBlLCByYXRlPXJhdGUpICogKDEtcCkpCgpzdGRpLnBsb3QgPC0gZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIikgJT4lIAogIHNlbGVjdChkYXksIHRlbXBzLCBjdW11bGF0aXZlX2dlcm0sIHRvdGFsX3NlZWRzKSAlPiUKICBtdXRhdGUodGVtcHM9YXMuZmFjdG9yKHRlbXBzKSwKICAgICAgICAgcHJvcF9nZXJtPWN1bXVsYXRpdmVfZ2VybS90b3RhbF9zZWVkcykKCnBvc3RfcGxvdDEuN212ICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9cHJvcF9nZXJtLGNvbG9yPXRlbXBzLGdyb3VwPXRlbXBzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGRhdGE9c3RkaS5wbG90KQogIApgYGAKCm1ha2UgYSBjb21wYXJpc29uIHBsb3QKCmBgYHtyfQpwb3N0X3Bsb3RfY29tcGFyZSA8LSBiaW5kX3Jvd3Mob2xkPXBvc3RfcGxvdDEuNywgbmV3PXBvc3RfcGxvdDEuN212LCAuaWQgPSAiZml0IikgCnBvc3RfcGxvdF9jb21wYXJlICU+JQogIGdncGxvdChhZXMoeD1kYXkseT1wcm9wX2dlcm0sY29sb3I9dGVtcHMsZ3JvdXA9c3RyX2ModGVtcHMsZml0KSxsaW5ldHlwZT1maXQpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoZGF0YT1zdGRpLnBsb3QsIGFlcyh4PWRheSx5PXByb3BfZ2VybSxjb2xvcj10ZW1wcyksIGluaGVyaXQuYWVzID0gRkFMU0UpCiAgCmBgYApvdmVyYWxsIG5vdCBhcyBnb29kCgojIyBleHRlbmQgdG8gbXVsdGlwbGUgc3BlY2llcy9wb3BzCgptb2RpZnkgbW9kZWwgMS43YQoKYGBge3J9Cmxlbmd0aCh1bmlxdWUoZ2VybW9uZSRwb3BzKSkgIzI0IHBvcHMKCmQgPC0gbGlzdChOPW5yb3coZ2VybW9uZSksCiAgICAgICAgICBnZXJtPWdlcm1vbmUkZ2VybSwKICAgICAgICAgIHRlbXBzPWFzLm51bWVyaWMoYXMuZmFjdG9yKGdlcm0uc3RkaSR0ZW1wcykpLAogICAgICAgICAgcG9wcz1hcy5udW1lcmljKGFzLmZhY3RvcihnZXJtb25lJHBvcHMpKSwKICAgICAgICAgIGRheT1nZXJtLnN0ZGkkZGF5KQoKc3Rhbm1vZGVsMi4xIDwtCiAgIgpkYXRhewogICAgaW50PGxvd2VyPTE+IE47ICAvLyBudW1iZXIgb2Ygb2JzZXJ2YXRpb25zCiAgICBpbnQgZ2VybVtOXTsKICAgIHZlY3RvcltOXSBkYXk7CiAgICBpbnQgdGVtcHNbTl07Cn0KcGFyYW1ldGVyc3sKICAgIHJlYWw8bG93ZXI9MD4gc2hhcGU7IC8vIHNob3VsZCBzZXQgbG93ZXIgYm91bmQKICAgIHZlY3Rvcls4XSBhOyAvLyBhbHBoYSBmb3IgdGhlIGdhbW1hIGN1cnZlLCBvbmUgZm9yIGVhY2ggdGVtcAogICAgdmVjdG9yWzhdIGFwOyAvLyBhbHBoYSBmb3IgdGhlIHByb3BvcnRpb24gZG9ybWFudCwgb25lIGZvciBlYWNoIHRlbXAKfQptb2RlbHsKICAgIHZlY3RvcltOXSBwOwogICAgdmVjdG9yW05dIG11OwogICAgYSB+IG5vcm1hbCggMCAsIC41ICk7IC8vbmFycm93IHByaW9ycyB0byBvdmVyY29tZSBkaXZlcmdlbnQgdHJhbnNpdGlvbnMKICAgIGFwIH4gbm9ybWFsKCAwLCAxLjUgKTsgCiAgICBzaGFwZSB+IGV4cG9uZW50aWFsKC41KTsgCiAgICBmb3IgKGkgaW4gMTpOKSB7CiAgICAgICAgcFtpXSA9IGFwW3RlbXBzW2ldXTsKICAgICAgICBwW2ldID0gaW52X2xvZ2l0KHBbaV0pOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSB7CiAgICAgICAgbXVbaV0gPSBhW3RlbXBzW2ldXTsKICAgICAgICAgLy8gYXBwbHkgdGhlIGludmVyc2UgbGluayBmdW5jdGlvbgogICAgICAgIG11W2ldID0gc2hhcGUgKiBleHAoLShtdVtpXSkpOwogICAgfQogICAgZm9yICggaSBpbiAxOk4gKSAKICAgICAgIGlmICggZ2VybVtpXSA9PSAwICkgdGFyZ2V0ICs9IGxvZ19taXgocFtpXSwgMCwgZ2FtbWFfbGNjZGYoZGF5W2ldIHwgc2hhcGUsIG11W2ldKSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSB0YXJnZXQgKz0gbG9nMW0ocFtpXSkgKyBnYW1tYV9scGRmKGRheVtpXSB8IHNoYXBlLCBtdVtpXSk7CiAgICBmb3IgKCBpIGluIDE6TiApIAogICAgICAgaWYgKCBnZXJtW2ldID09IDEgKSBkYXlbaV0gfiBnYW1tYSggc2hhcGUsIG11W2ldICk7Cn0KIgptMS43YSA8LSBzdGFuKG1vZGVsX2NvZGU9c3Rhbm1vZGVsMS43YSwgZGF0YT1kLCBjaGFpbnM9NCwgY29yZXM9NCwgY29udHJvbD1saXN0KGFkYXB0X2RlbHRhPS45OSkpCmBgYAoKCgoKIyMgaG93IGNhbiB3ZSBzd2l0Y2ggdGhpcyBhcm91bmQgdG8gdXNlIGRheSBhcyBhIHByZWRpY3Rvcj8K